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Computational Chemistry 
Glossary 



• Ab initio. The latin term "ab initio" means "from the beginning" or "from 
first principles". Within ab-initio simulations the electrons are treated quan- 
tum mechanically. If we are interested in describing electron distribution in 
detail there is no substitute for quantum mechanics. If we don't make use of 
any (experimentally derived) parameter in solving the Schrodinger equation 
for the electrons, the theory is called ab initio. If the nuclei are described 
at classical level and the electrons at a quantum mechanical level the ap- 
proach is called semiclassical. In this course we will always remain in this 
approximation. 

• Approximations in computational quantum chemistry 

We distinguish two main types of approximations: 

- either in the Hamiltonian (for instance by changing from a a wavefunc- 
tion based to a density based description of the electronic interaction 
(DFT), or through the simplification of the electronic interaction term 
in semi-empirical, wavefunction based models). 

- or in the description of the many-electron wavefunction. 

In computational quantum chemistry the electronic wavefunction of a molec- 
ular system is often expanded as a sum of anti-symmetrized many electron 
wavefunctions (Slater determinants) 

^ e i(r 1 ,si,r 2 ,. ■ ■ ,r N ,s N ) = 



The components of the Slater determinant, mj (rj), are one-electron molec- 
ular orbitals which are usually given as an expansion in "atomic orbitals", 
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(here r stays for the Cartesian coordinates (x, y, z) and s is the spin of the 
electron (s € {a,f3}). 

The collection of coefficients D ., and C... fully characterizes the solution of 
the electronic Schrodinger equation for atoms and molecules. 

The most commonly used approximate methods for the solution of the elec- 
tronic molecular Schrodinger equation are: 

- Semi-empirical (MNDO, AMI, PM3, etc.): use a single Slater determi- 
nant (only one C is 1 all the others are equal 0). Vary the coefficients D, but 
just use empirical estimates rather than the true integrals. 

Very cheap, but only accurate for molecule similar to those used to develop 
the empirical estimates. 

- DFT (B3LYP, BLYP, PW91, etc.): slightly empirical, but much more reli- 
able than semi-empirical methods. CPU: cheap, same as HF 0(A 3 ). Errors 
~ 4 kcal/mole (comparable accuracy to MP2 but much cheaper). Preferred 
method for geometries, second derivatives and transition-metal containing 
systems. 

- HF (Hartree-Fock, SCF): only one many-electrons Slater determinant is 
used. Vary the D's, all terms in the electronic Hamiltonian calculated 'ab- 
initio' within the mean field approximation, no empirical parameters. 
CPU: cheap 0(N 3 ) errors ~ 15 kcal/mol. 

- MP2, MP4 (Moller-Plesset, MBPT): Vary the D's first, then set the C's 
to the values given by perturbation theory (you don't freely vary these C's, 
saving CPU). 

MP2: medium CPU: 0(N 5 ), errors ~ 5 kcal/mol. 

- CI, CISD, QCISD (Configuration Interaction): Vary the coefficients D 
first, freeze them, then vary a lot of the coefficients C. 

Expensive. Not used much anymore, CCSD is preferred. 

- MCSCF, CASSCF: vary a finite set of C's and all the D's simultaneously. 

Expensive. Good for understanding cases where several electronic states have 
comparable energies. User expertise required to select which C's to vary. 

- CAS-PT2: Determine the D's and some C's by CASSCF, then determine 
more C's by perturbation theory. 
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Not much more expensive than CASSCF. Sometimes very good, but not reli- 
able. 

- MRCI (multi reference CI): Determine the D's and some C's by CASSCF 
or MCSCF, freeze these, then allow many of the C's to vary. 

Super expensive. Very high accuracy for small systems. 

- CCSD, CCSD(T) (Coupled Cluster): Vary the D's, fix them, then vary 
a lot of the C's, but constraining certain relationships between the C's. This 
allows you to effectively use a longer expansion without increasing the num- 
ber of adjustable parameters so much. The constraints force the solution 
to be "size-consistent", i.e. two molecules calculated simultaneously have 
exactly the same energy as two molecules calculated separately. 
Expensive. Often very accurate. 

- Extrapolations ("Composite Methods"): G2, G3, CBS-q, CBS-Q, CBS- 
QB3, CBS-RAD Run a series of the above calculations with different size 
basis sets, following some recipe. The results from all these calculations are 
extrapolated to an estimate of the true potential V{R). These methods give 
excellent accuracy in less CPU time than CCSD or MRCI. However, the mul- 
tiple steps involved provide many opportunities for something to go wrong. 
Accuracy: usually 1-2 kcal/mol. 

• Ab initio forces. Among the four fundamental forces in nature (electro- 
static, weak and strong interactions, and gravitation) the only one that is 
relevant to ab initio simulations is the Coulomb interaction among nuclei 
and electrons. 

• Atom. In chemistry and physics, an atom (Greek for "uncuttable") is the 
smallest possible particle of a chemical element that retains its chemical 
properties. Whereas the word atom originally denoted a particle that cannot 
be cut into smaller particles, the atoms of modern parlance are composed of 
subatomic particles: 

— electrons, which have a negative charge and are the least massive of the 
three; 

— protons, which have a positive charge and are about 1836 times more 
massive than electrons; and 

— neutrons, which have no charge and are about 1839 times more massive 
than electrons. 

(Protons and neutrons are not fundamental particels, but they are composed 
of quarks bound together by gluons (strong interaction, confinement.) Atoms 
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can differ in the number of each of the subatomic particles they contain. 
Atoms of the same element have the same number of protons (called the 
atomic number). Within a single element, the number of neutrons may vary, 
determining the isotope of that element. The number of electrons associated 
with an atom is most easily changed, due to the lower energy of binding of 
electrons. The number of protons (and neutrons) in the atomic nucleus may 
also change, via nuclear fusion, nuclear fission or radioactive decay, in which 
case the atom is no longer the same element it was. This requires much 
larger energies than changes in the electronic composition. 
Atoms are electrically neutral if they have an equal number of protons and 
electrons. Atoms which have either a deficit or a surplus of electrons are 
called ions. Electrons that are furthest from the nucleus may be transferred 
to other nearby atoms or shared between atoms. By this mechanism atoms 
are able to bind into molecules and other types of chemical compounds like 
ionic and covalent network crystals. 

In all the methods we discuss in this course, the nucleus is considered as 
a classical point charge with no extension, carrying a charge Zj. Its position 
in space is defined by the position vector Rj. The electrons are described 
quantum mechanically. Sometimes, core electrons (the inner electrons in 
the shells close to the nucleus) are not considered explicitly because they 
are mainly chemically inert. In this case, the nuclear charge is replaced by, 
Zj — N c , where N c is the number of core electrons, and the effects of the core 
electrons on the remaining valence electrons is replaced by the corresponding 
pseudopotentials. 

• Classical mechanics. Classical mechanics is concerned with the set of phys- 
ical laws governing and mathematically describing the motions of bodies and 
aggregates of bodies in the classical limit (neglecting quantum effects) . The 
term classical mechanics was coined in the early 20th century to describe the 
system of mathematical physics developed in the 400 years since the ground- 
breaking works of Brahe, Kepler, and Galileo, but before the development 
of quantum physics and relativity. The initial stage in the development of 
classical mechanics is often referred to as Newtonian mechanics, and is associ- 
ated with the mathematical methods invented by Newton himself, in parallel 
with Leibniz, and others. More abstract, and general methods include La- 
grangian mechanics and Hamiltonian mechanics. While the terms classical 
mechanics and Newtonian mechanics are usually considered equivalent, the 
conventional content of classical mechanics was created in the 19th century 
and differs considerably (particularly in its use of analytical mathematics) 
from the work of Newton. 

• Computational chemistry. Computational chemistry is a branch of chem- 
istry that uses the results of theoretical chemistry incorporated into efficient 
computer programs to calculate the structures and properties of molecules 



10 



and solids, applying these programs to real chemical problems. Examples 
of such properties are molecular structures (i.e. the expected positions of 
the constituent atoms), energy and interaction energy, charges, dipoles and 
higher multipole moments, vibrational frequencies, reactivity and spectro- 
scopic quantities. The term computational chemistry is also sometimes used 
to cover any of the areas of science that overlap between computer science 
and chemistry. Electronic structure theory is a subdiscipline of computa- 
tional chemistry. 

n-body and many-body problems. The "n-body problem" is the prob- 
lem of finding, given the initial positions, masses, and velocities of n bodies, 
their subsequent motions as determined by classical mechanics, i.e., Newton's 
laws of motion. The equivalent problem in quantum mechanics is named 
"many-body problem". 

In classical mechanics (n-body problem) an exact solution is only known for 
the case n = 2. The three-body problem (n = 3) is much more complicated; 
its solution can be chaotic. 

In quantum mechanics the "many-body problem" is usually posed as the 
question of solving more complex problems than the hydrogen atom. De- 
pending on the complexity of the molecule, different models are used. The 
kind of methods and approximations that can be used is the main 
subject of this course. 

Force field. In Force Field Methods the energy of a molecular system is 
written as a parametric function of the nuclear coordinates. The parameters 
that enter the function are fitted to experimental or higher level computa- 
tional data. The molecules are modelled as atoms held together by bonds. 
Three and four atom energy terms are also included to describe angles and 
dihedrals. All atoms are represented by point charges with fractional charges 
that interact through a Coulomb potential. The ensemble of all parameters 
is called a Force Field. For a given molecule, the assignment of force field 
parameters to the differend atoms, bonds, angles, and dihedrals defines its 
topology. Once the topology is defined the constraints to the dynamics of the 
system are set once forever. No chemical reactions can occur. 
Forces on the atoms are computed as derivative of the total energy and the 
system evolves according to Newton's law. 

Fundamental forces. Traditionally, modern physics describes all natural 
phenomena using four fundamental interactions: gravitation, electromag- 
netism, the weak interaction, and the strong interaction. 
In this course we consider exclusively electrostatic interactions between 
nuclei and electrons. 

Hamiltonian operator. In quantum mechanics, the quantum Hamiltonian 
H is the operator that determines to the total energy of the system, i.e. the 
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total energy is the observable associated with H. As with all observables, the 
spectrum of the Hamiltonian are the possible outcomes when one measures 
the total energy of a system. 

• Molecular dynamics (MD) and ab-initio molecular dynamics. Once 
the Hamiltonian is given (classical or semiclassical) one can compute its 
derivative with respect to the nuclear positions to calculate the forces that 
act on the nuclei. These are used together with the Newton equation to 
compute the dynamics of the nuclei. Nuclear dynamics can be carried out in 
different thermodynamic ensembles, microcanonical (NVE, constant energy), 
canonical (NVT or NPT, constant temperature or grand canonical (/zVT, 
variable number of particles). 

The electronic dynamics can be carried out explicitly by solving the time- 
dependent Schrodinger equation, or - due to the much smaller mass of the 
electrons - can be considered as an instantaneous relaxation of the electrons 
into the new energy minimum provided by the new atomic positions. In the 
second case, named Born-Oppenheimer dynamics, there is no real electron 
dynamics , but the only dynamical variables are the nuclear positions. 

MjRjit) = V / min{($ |H( J R/)|$o)} 
$o 

Born-Oppenheimer MD can be accelerated using the Car-Parrinello MD 
scheme. 

• Molecular mechanics. The term molecular mechanics refers to the use 
of Newtonian mechanics to model molecular systems. Molecular mechanics 
approaches are widely applied in molecular structure refinement, molecular 
dynamics simulations, Monte Carlo simulations and ligand docking simula- 
tions. Molecular mechanics can be used to study small molecules as well 
as large biological systems or material assemblies with many thousands to 
millions of atoms. All-atomistic molecular mechanics methods have the fol- 
lowing properties: 

— Each atom is simulated as a single hard spherical particle. 

— To each such particle a radius (typically the van der Waals radius) 
and a constant net charge (generally derived from high-level quantum 
calculations and/or experiment) is assigned. 

— Bonded interactions are treated as "springs" with equilibrium distances 
and harmonic force constants equal to the experimental or calculated 
bond lengths and vibrational frequencies. 

• Optimization of wavefunctions and geometries. By wavefunction op- 
timization we refer to the proccess of finding numerical approximations of 
the system's wavefunction by successive iterations. This process requires 
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a variational principle that guarantees convergence by minimization of the 
total energy. 

E < . 

During a geometry optimization, nuclear forces are computed at the end 
of each wavefunction optimization process. The nuclei are then shifted in 
direction of the computed forces and a new wavefunction is computed for 
the new positions. The process is iterated until convergence. The final 
nuclear coordinates correspond to the (global) minimum of the potential 
energy surface. 

Pseudopotentials. In quantum mechanics, the pseudopotential formal- 
ism is an attempt to replace the complicated effects of the motion of the 
core (i.e. non-valence) electrons of an atom or ion and its nucleus with an 
effective potential, or pseudopotential, so that the Schrodinger equation con- 
tains a modified potential term instead of e.g. the Coulombic potential term 
normally found in the Schrodinger equation and the electronic Schrodinger 
equation can be solved for the valence electrons only. 

Quantum chemistry. Quantum chemistry is a branch of theoretical chem- 
istry, which applies quantum mechanics to address issues and problems in 
chemistry. The description of the electronic behavior of atoms and molecules 
to describe their reactivity is one of the applications of quantum chemistry. 
Quantum chemistry lies on the border between chemistry and physics, and 
significant contributions have been made by scientists from both fields. It has 
a strong and active overlap with the fields of atomic physics and molecular 
physics, as well as physical chemistry. 

Quantum mechanics. Quantum mechanics is a fundamental branch of 
theoretical physics that replaces classical mechanics and classical electromag- 
netism at the atomic and subatomic levels. It is the underlying mathemati- 
cal framework of many fields of physics and chemistry, including condensed 
matter physics, atomic physics, molecular physics, computational chemistry, 
quantum chemistry, particle physics, and nuclear physics. Along with general 
relativity, quantum mechanics is one of the pillars of modern physics. 

Relativistic quantum mechanics. Relativistic effects in quantum me- 
chanics are described by the Dirac equations. The notion of wavefunction 
and spin is replaced by the four-components spinor (solution of the Dirac 
equation). In the non-relativistic limit 

e~^« P c (2) 

the Dirac equations reduces to the Schrodinger equation. Here p is the 
particle momentum, m is its mass and c is the speed of light. 
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In this course we do not consider relativistic effects. 

• Schrodinger equation. In physics, the Schrodinger equation, proposed by 
the Austrian physicist Erwin Schrodinger in 1925, describes the space- and 
time-dependence of quantum mechanical systems. It is of central importance 
to the theory of quantum mechanics, playing a role analogous to Newton's 
law in classical mechanics. 

In the mathematical formulation of quantum mechanics, each system is as- 
sociated with a complex Hilbert space such that each instantaneous state 
of the system is described by a unit vector in that space. This state vector 
encodes the probabilities for the outcomes of all possible measurements ap- 
plied to the system. As the state of a system generally changes over time, 
the state vector is a function of time. The Schrodinger equation provides a 
quantitative description of the rate of change of the state vector. 

• Statistical Mechanics. Statistical mechanics is the application of proba- 
bility theory, which includes mathematical tools for dealing with large pop- 
ulations, to the field of mechanics, which is concerned with the motion of 
particles or objects when subjected to a force. It provides a framework for 
relating the microscopic properties of individual atoms and molecules to the 
macroscopic or bulk properties of materials that can be observed in everyday 
life, therefore explaining thermodynamics as a natural result of statistics and 
mechanics (classical and quantum) at the microscopic level. In particular, 
it can be used to calculate the thermodynamic properties of bulk materials 
from the spectroscopic data of individual molecules. 
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Chapter 1 



Introduction to computational 
quantum chemistry 

Computational quantum chemistry is a branch of theoretical chemistry whose 
major goal is to create efficient mathematical approximations and computer 
programs to calculate the properties of molecules such as their total energy, 
dipole and quadrupole moments, potential energy (free energy) sur- 
faces, vibrational frequencies, excitation energies and other diverse 
spectroscopic quantities, reactivity behavior, the involved reaction 
mechanism and reaction dynamics. The term is also sometimes used to 
cover the areas of overlap between computer science and chemistry. 

Theoretical vs. computational chemistry. The term theoretical chem- 
istry may be defined as a mathematical description of chemistry, whereas 
computational chemistry is usually used when a mathematical method is 
sufficiently well developed that it can be automated for implementation on 
a computer. Note that the word exact does not appear here, as very few 
aspects of chemistry can be computed exactly. Almost every aspect of chem- 
istry, however, can be and has been described in a qualitative or approximate 
quantitative computational scheme. 

Accuracy vs. efficiency. It is, in principle, possible to use one very accu- 
rate method and apply it to all molecules. Although such methods are well- 
known and available in many computer programs, the computational cost 
of their use grows exponentially with the number of electrons. Therefore, 
a great number of approximate methods strive to achieve the best trade-off 
between accuracy and computational cost. Present computational chemistry 
can routinely and very accurately calculate the properties of molecules that 
contain no more than 40-50 electrons. The treatment of molecules that con- 
tain up to l'OOO electrons is computationally tractable by DFT. 

The main target of this course is to introduce the basic notions 
of computational quantum chemistry, which allow to explore reac- 
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tion mechanisms and explain observations of laboratory reactions. 

Several major areas may be distinguished to which computational quan- 
tum chemistry methods can be applied: 

• The computation of physical and chemical properties of molecules (struc- 
ture and dynamics) 

• Identifying correlations between chemical structures and properties 

• Understanding reaction mechanisms and thermodynamic properties 

• Computation of the potential energy surfaces. Determination of quan- 
tum forces on the nuclei and correspondent ab initio molecular dynam- 
ics. 

• Computational approaches to help in the efficient synthesis of com- 
pounds 

• Computational approaches to design molecules that interact in specific 
ways with other molecules (e.g. drug design) 

Methods for quantum chemistry can also be applied to solid state physics 
problems. The electronic structure of a crystal is in general described by a 
band structure, which defines the energies of electron orbitals for each point 
in the Brillouin zone. 



1.1 Ab Initio Methods 

The programs used in quantum chemistry are based on different quantum- 
mechanical methods that solve the molecular Schrodinger equation associated 
with the molecular Hamiltonian. Methods that do not include empirical or 
semi-empirical parameters in their equations - are derived directly from theo- 
retical principles, with no inclusion of experimental data - are generally called 
ab initio methods. However, the same term is also used to design theories 
that are derived from exact quantum mechanical equations but which also 
assume a certain level of approximation. The approximations made in these 
cases are usually mathematical in nature, such as using a simpler functional 
form or getting an approximate solution for a complicated differential equa- 
tion. 



1.1.1 Electronic structure 

One of the primary goals of quantum chemical methods is to determine the 
electronic structure, e.g. the probability distribution of electrons in chemical 
systems. The electronic structure is determined by solving the Schrodinger 
equation associated with the electronic molecular Hamiltonian. In this pro- 
cess, the molecular geometry is considered as a fixed parameter. Once the 
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optimal electronic wavefunction is determined, one can access the gradient 
on each nuclei (as the derivative of the total energy with respect to the posi- 
tions of the nuclei) and update their positions accordingly until the process 
reaches convergence. This is how we obtain optimized molecular geometries. 
Usually the basis set (which is usually built from the LCAO ansatz) used to 
solve the Schrodinger equation is not complete and does not span the full 
Hilbert space associated to the system. However, this approximation allows 
one to treat the Schrodinger equation as a "simple" eigenvalue equation of 
the electronic molecular Hamiltonian with a discrete set of solutions. The 
problem of dealing with function (functionals) and operators can be there- 
fore translated into a linear algebra calculation based on energy matrices 
and state vectors. The obtained eigenvalues are functions of the molecular 
geometry and describe the potential energy surfaces. 
Many optimized linear algebra packages have been developed for this purpose 
(e.g., LAPACK - Linear Algebra PACKage: http:/ /www. netlib.org/lapack). 

The most common type of an ab initio electronic structure approach is 
called the Hartree-Fock (HF) method, in which the Coulombic electron- 
electron repulsion is taken into account in an averaged way (mean field ap- 
proximation). This is a variational calculation, therefore the obtained ap- 
proximate energies, expressed in terms of the system's wavefunction, are 
always equal to or greater than the exact energy, and tend to a limiting 
value called the Hartree-Fock limit. Many types of calculations begin with 
a HF calculation and subsequently correct for the missing electronic correla- 
tion. M0ller-Plesset perturbation theory (MP) and Coupled Cluster 
(CC) are examples of such methods. 

A method that avoids making the variational overestimation of HF in the 
first place is Quantum Monte Carlo (QMC), in its variational, diffusion, 
and Green's functions flavors. These methods work with an explicitly corre- 
lated wavefunction and evaluate integrals numerically using a Monte Carlo 
integration. Such calculations can be very time consuming, but they are 
probably the most accurate methods known today. 

Density Functional Theory (DFT) methods are often considered to 
be ab initio methods for determining the molecular electronic structure, even 
though one utilizes functionals usually derived from empirical data, proper- 
ties of the electron gas, rare gas calculations or more complex high level 
approaches. In DFT, the total energy is expressed in terms of the total 
electron density, rather than the wavefunction. 

The Kohn-Sham formulation of DFT introduces a set of non-interacting 
fictitious molecular orbitals, which makes the chemical interpretation of DFT 
simpler, even though their nature is different from the rea/-HF ones. In addi- 
tion, KS based DFT has the advantage to be easily translated into an efficient 
computer code. 

The most popular classes of ab initio electronic structure methods are: 
• Hartree-Fock 
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• Generalized Valence Bond 

• M0ller-Plesset perturbation theory 

• Multi-Configurations Self Consistent Field (MCSCF) 

• Configuration interaction 

• Multi-Reference Configuration Interaction 

• Coupled cluster 

• Quantum Monte Carlo 

• Reduced density matrice approaches 

• Density functional theory 

Ab initio electronic structure methods have the advantage that they can 
be made to converge systematically to the exact solution (by improving the 
level of approximation) . The convergence, however, is usually not monotonic, 
and sometimes - for a given property - a less accurate calculation can give 
better result (for instance the total energy in the series HF, MP2, MP4, . . . ). 

The drawback of ab initio methods is their cost. They often take enor- 
mous amounts of computer time, memory, and disk space. The HF method 
scales theoretically as N A (N being the number of basis functions) and cor- 
related calculations often scale much less favorably (correlated DFT calcula- 
tions being the most efficient of this lot). 



1.1.2 Chemical dynamics 

Once the electronic and nuclear motions are separated (within the Born- 
Oppenheimer approximation), the wave packet corresponding to the nuclear 
degrees of freedom can be propagated via the time evolution operator asso- 
ciated with the time- dependent Schrodinger equation (for the full molecular 
Hamiltonian) . The most popular methods for propagating the wave packet 
associated to the molecular geometry are 

• the split operator technique for the time evolution of a phase-space 
distribution of particles (nuclei) subject to the potential generated by 
the electrons (Liouville dynamics), 

• the Multi-Configuration Time-Dependent Hartree method (MCTDH), 
which deals with the propagation of nuclear-wavepackets of molecular 
system evolving on one or several coupled electronic potential energy 
surfaces. 

• Feynman path-integral methods. In this approach the quantum par- 
tition function describing the dynamics of N nuclei is mapped into a 
classical configurational partition function for a N x P-particle system, 
where P are discrete points along a cyclic path. 
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The semiclassical approach. In this case one adopts a classical molec- 
ular mechanics propagation of the nuclear dynamics (where the equation 
of motion is given by Newton's law) combined with any treatment of the 
time- dependent Schrodinger equation for the electrons. The combination of 
classical molecular dynamics with an efficient scheme for the computation of 
the potential energy surface, like for instance DFT, allows for ab inito MD 
simulations of systems of the size of few hundreds of atoms. Non-adiabatic 
effects 1 can also be introduced by means of the fewest switches surface hop- 
ping scheme of Tully et. al. or the mean-field approximation (Ehrenfest 
dynamics) . 

During this course we will only deal with the semiclassical approximation. 

Different numerical schemes for the update of the electronic state at each 
molecular step are possible 

• uncorrelated wavefunction optimization at the new nuclear configura- 
tion. This method is often called Born-Oppenheimer Molecular Dy- 
namics (BO MD). 

• within DFT, the electronic density can be propagated using the "ficti- 
tious" electronic dynamics obtained from the Car-Parrinello extended 
Lagrangian scheme. In this way one saves the cost a full wavefunc- 
tions optimization at each MD step. The method is referred as Car 
Parrinello Molecular Dynamics (CP MD). 

• propagation of the electronic wavefunction using an approximate so- 
lution of the time-dependent Schrodinger equation. Due to the small 
mass of the electrons (compared to the one of the nuclei) the time step 
for the propagation of the electrons is necessary much smaller than the 
one required for the solution of the Newton equation for the ions. 



1.2 Semiempirical methods 
1.2.1 Electronic structure 

The most expensive term in a Hartree-Fock calculation, are the so-called 
two-electron integrals. Semiempirical methods are based on the HF method 
but the costly two-electrons integrals are approximated by empirical data or 
omitted completely. In order to correct for this approximation, semiempirical 
methods introduce additional empirical terms in the Hamiltonian which are 
weighted by a set of a priori undefined parameters. In a second step, these 
are fitted in order to reproduce results in best agreement with experimental 
data. 



1 Non-adiabatic effects describe the quantum phenomena that arise when the dynamics 
is approaching regions in phase space where two or more adiabatic PES are approaching 
together. 
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Semiempirical calculations are much faster than their ab initio counter- 
parts. Their results, however, can be badly wrong if the molecule being 
computed is not similar enough to the molecules in the database used to 
parameterize the method. 

Semiempirical calculations have been very successful in the description 
of organic chemistry, where only a few elements are used extensively and 
molecules are of moderate size. 

Semiempirical methods also exist for the calculation of electronically ex- 
cited states of polyenes. These methods, such as the Pariser-Parr-Pople 
(PPP) method, can provide good estimates of the electronic excited states, 
when parameterized well. Indeed, for many years, the PPP method outper- 
formed ab initio excited state calculations. Large scale ab initio calculations 
have confirmed many of the approximations of the PPP model and explain 
why the PPP-like models work so well in spite of its simple formulation. 

1.3 Mixed Quantum Mechanical/Molecular me- 
chanical (QM/MM) Methods 

The combination of ab initio methods (QM) with a classical molecular me- 
chanical (MM) treatment of the environment enables the investigation of 
larger systems in which the chemical active region (QM) is confined in a 
localized volume surrounded by MM atoms. This approach is particularly 
suited for the analysis of enzymatic reactions in proteins. 



1.4 Software packages 

A number of software packages that are self-sufficient and include many 
quantum-chemical methods are available. The following is a table (not com- 
plete) illustrating the capabilities of various software packages: 



Package 


MM 


Semi- Empirical 


HF 


Post-HF 


DFT 


Ab-inito MD 


Periodic 


QM 'MM 


ACES 


N 


N 


Y 


Y 


N 


N 


N 


N 


ADF 


N 


N 


N 


N 


Y 


N 


Y 


Y 


CPMD 


Y 


N 


N 


N 


Y 


Y 


Y 


Y 


DALTON 


N 


N 


Y 


Y 


Y 


N 


N 


N 


GAUSSIAN 


Y 


Y 


Y 


Y 


Y 


Y(?) 


Y 


Y 


GAMESS 


N 


Y 


Y 


Y 


Y 


N 


N 


Y 


MOLCAS 


N 


N 


Y 


Y 


N 


N 


N 


N 


MOLPRO 


N 


N 


Y 


Y 


Y 


N 


N 


N 


MOPAC 


N 


Y 


N 


N 


N 


N 


Y 


N 


NWChem 


Y 


N 


Y 


Y 


Y 


Y(?) 


Y 


N 


PLATO 


Y 


N 


N 


N 


Y 


N 


Y 


N 


PSI 


N 


N 


Y 


Y 


N 


N 


N 


N 


Q-Chem 


? 


N 


Y 


Y 


Y 


N 


N 


N 


TURBOMOL 


N 


N 


Y 


Y 


Y 


Y 


Y 


N 
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Chapter 2 



Basic Principles of quantum 
mechanics 

For a repetition of quantum mechanics you may find the following links use- 
ful: 

http:/ / en.wikipedia. org/wiki/ Quantum _ mechanics 
http:/ / ocw. mil. edu/ Ocw Web / Chemistry/index, htm 
http://ocw.mit. edu/ Ocw Web /Physics/index, htm 
http:/ /bethe. Cornell, edu/ 
http:/ /plato. Stanford, edu/ entries / qm 



2.1 Postulates of Quantum Mechanics 

• Postulate 1. The state of a quantum mechanical system is completely 
specified by a function ty(r,t) that depends on the coordinates of all 
particles and on time. This function, called the wave function or state 
function, has the important property that \l/*(r, (r, t)dr is the prob- 
ability that the particle lies in the volume element dr located at r at 
time t. The wavefunction must satisfy certain mathematical conditions 
because of this probabilistic interpretation. For the case of a single 
particle, the probability of finding it somewhere in space is 1, so that 
we have the normalization condition 

/oo 
ty*(r,t)^(r,t)dT = 1 (2.1) 
-oo 

It is customary to also normalize many-particle wavefunctions to 1. 
The wavefunction must also be single- valued, continuous, and finite. 

• Postulate 2. To every observable in classical mechanics there corre- 
sponds a linear, Hermitian operator in quantum mechanics. If we re- 
quire that the expectation value of an operator A is real, then A must 
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be a Hermitian operator. Some common operators occurring in quan- 
tum mechanics are collected in the following table. 



Observable 


Observable 


Operator 


Operator 


Name 


Symbol 


Symbol 


Operation 


Positron 


_r 


r 


Multiply by JL 


Momentum 


P 


P 




Kinetic energy 


T 


T 




Potential energy 




Ha 


Multiply by V{t) 


Total energy 


E 


k 




Angular 
momentum 


k 


t, 






k 










I 





Figure 2.1: Physical observables and their corresponding 
quantum operators. 



• Postulate 3. In any measurement of the observable associated with op- 
erator A, the only values that will ever be observed are the eigenvalues 
a, which satisfy the eigenvalue equation 

AV(r,t) = a*(r,i) (2.2) 

This postulate captures the central point of quantum mechanics-the 
values of dynamical variables can be quantized (although it is still pos- 
sible to have a continuum of eigenvalues in the case of unbound states). 

If the system is in an eigenstate of A with a single eigenvalue a, then 
any measurement of the quantity A will yield a. 

Although measurements must always yield an eigenvalue, the state does 
not have to be an eigenstate of A initially. An arbitrary state can 
be expanded in the complete set of eigenvectors of A (A^i(r,t) = 
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ai^i(r,t)) as 



N 



(2.3) 



i=X 



where N may go to infinity. In this case we only know that the mea- 
surement of A will yield one of the values Oj, but we don't know which 
one. However, we do know the probability that eigenvalue a, will occur: 
it is the absolute value squared of the coefficient, |cj| 2 , leading to the 
fourth postulate below. 

An important second half of the third postulate is that, after measure- 
ment of \l/ yields some eigenvalue a^, the wavefunction immediately 
"collapses" into the corresponding eigenstate ^j. Thus, measurement 
affects the state of the system. This fact is used in many elaborate 
experimental tests of quantum mechanics. 



• Postulate 4. If a system is in a state described by a normalized wave- 
function then the expectation value of the observable corresponding 

to A is given by 



• Postulate 5. The wavefunction or state function of a system evolves in 
time according to the time- dependent Schrodinger equation 



The central equation of quantum mechanics must be accepted as a 
postulate. 



• Postulate 6. The total wavefunction must be antisymmetric with re- 
spect to the interchange of all coordinates of one fermion 1 with those 
of another. Electronic spin must be included in this set of coordinates. 
The Pauli exclusion principle is a direct result of this antisymmetry 
principle. We will later see that Slater determinants provide a conve- 
nient means of enforcing this property on electronic wavef unctions. 



1 Fermins: particles with half-integer spins. Electrons are fermins. 
Bosons: particles with integer spin, e.g. the nucleous of a C-12 atom. 




(2.4) 




(2.5) 
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2.2 The molecular Hamiltonian and the Born- 
Oppenheimer approximation 

This chapter was adapted from the lecture notes "A Short Summary of Quan- 
tum Chemistry", MITOPENCOURSEWARE, December 2004. 

Quantum Chemistry is (typically) based on the non-relativistic Schrodinger 
equation within the Born-Oppenheimer approximation. The Schrodinger 
equation is (we use atomic units: h = 1, m e ; ec = 1, e e \ ec = 1) 

H tot (R, P, r, p)*tot(R, P, r, p) = E(R, P) V tot (R, P, r, p) (2.6) 

where r, p = d/dr are the electronic collective coordinates and R, P = d/dH 
are the nuclear collective coordinates, and 

- E is an allowed energy of the system (the system is usually a molecule) . 

- ^>tot is a function of the positions of all the electrons and nuclei (we 
drop all spin dependencies). 

- H tot is a differential operator constructed from the classical Hamiltonian 
H(P, R, p, r) = E by replacing all the momenta p^ (resp.Pj) with 
{i)d/dvi ((i)d/dHi) as long as all the p (P) and r (R) are Cartesian. 

For a system of nuclei and electrons in vacuum with no external fields, 
neglecting magnetic interactions, using atomic units: 

/ n I<J In n<m 

(2.7) 

The Born-Oppenheimer approximation is to neglect some of the terms cou- 
pling the electrons and nuclei, so one can write: 

^(R, r) = V nucl (R) ^ eZec (r; R) (2.8) 

and 

Htot = f nud (P, R) + H elec (p, r; R) (2.9) 

which ignores the dependence of H e i ec on the momenta of the nuclei P. One 
can then solve the Schrodinger equation for the electrons (with the nuclei 
fixed, indicated by (; R)). The energy we compute will depend on the posi- 
tions R of those fixed nuclei, call it E(R): 

H elec {p, r; R)^ eiec (r; R) = E(R) ^ eZec (r; R) (2.10) 

The collection of all possible nuclear configurations, R together with the 
associated energies, E(R), defines a potential energy surface, V^(R) for the 
nuclei. 
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Now we can go back to the total Hamiltonian, and integrate over all the elec- 
tron positions r, ignoring any inconvenient term, to obtain an approximate 
Schrodinger equation for the nuclei: 



{V elec (r,R)\H tot \q elec (r,R)) = H nud = f nucl (P,R) + V(R) (2.11) 



2.2.1 The nuclear Schrodinger equation 

Both approximate Schrodinger equations (for electrons eq. 2.10 and for nuclei 
eq. 2.12) are still much too hard to solve exactly (they are partial differential 
equations in 3N particle coordinates), so we have to make more approxima- 
tions. 

V(R) is usually expanded to second order R about a stationary point Rq: 



and then the translations, rotations, and vibrations are each treated sepa- 
rately, neglecting any inconvenient terms that couple the different coordi- 
nates. In this famous "rigid-rotor-harmonic-oscillator (RRHO)" approxima- 
tion, analytical formulas are known for the energy eigenvalues, and for the 
corresponding partition functions Q (look in any Phys.Chem. textbook). 
This approximate approach has the important advantage that we do not 
need to solve the Schrodinger equation for the electrons at very many R's: 
we just need to find a stationary point R , and compute the energy and the 
second derivatives at that R . Many computer programs have been written 
that allow one to compute the first and second derivatives of V^(R) almost 
as quickly as you can compute V. For example, for the biggest calculation 
called for in problem 3 with 10 atoms and 3 * 10 = 30 coordinates Rj, it 
takes about half a minute on a PC to compute V(Ro) and only about 13 

more minutes to compute the 30 * 30 = 900 second derivatives ( |r~|1~ ) ■ If 



you tried to do this naively by finite differences, it would take about 15 hours 
to arrive at the same result (and it would probably be less accurate because 
of finite differencing numerical errors.) The analytical first derivatives are 
used to speed the search for the stationary point (e.g. the equilibrium ge- 
ometry) R . Often the geometry and the second derivatives are calculated 
using certain approximations, but the final energy V(Ro) is computed more 
accurately (since thermodynamics data and reaction rates are most sensitive 
to errors in V(Ro), and even poor approximations often get geometry and 
frequencies close to correct). 

Therefore, as long as a second-order Taylor expansion approximation for V 
is adequate we are in pretty good shape. Molecules and transition states 
with "large amplitude motions" (i.e. the Taylor expansion is not adequate) 



with 
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are much more problematic, dealing with them is an active research area. 
Fortunately, there are many systems where the conventional second-order V, 
RRHO approximation is accurate. 



2.2.2 The electronic Schrodinger equation 

The question is now how to compute the required potential V(R) which acts 
on the nuclei at a given geometry R. What we need to solve is 2.10, 

H dec (p, r; R)^ eZec (r; R) = V(R)^ dec (r; R) 
where in a vacuum, in the absence of fields, and neglecting magnetic effects 



H, 



elec' 



^E^+E^-E^+E^ 

n I<J 1 1 J 1 In 1 1 n ' n<m 1 m 



(2.14) 

and because the electrons are indistinguishable fermions any permutation of 
two electrons must change the sign of the wavefunction \E^ erf (r; R) (this is 
a really important constraint called the Pauli exclusion principle, it is the 
reason for the specific structure of the periodic table). 

In addition, because the spin is a good quantum number we can chose the 
electronic wavefunction to be simultaneously an eigenfunction of the spin 
operator: 



S 2 \V elec ) = S(S + l)\y elec ) 

$z I ^ dec) 'S'z | ^eiec) 



(2.15) 
(2.16) 



We can write ^ e iec in a form that will guarantee it satisfies the Pauli principle, 
namely using the Slater determinant many-electron wavef unctions: 



* ei (ri,r 2 , ...,r N ) 
where 

(r 2 ) ...(p. 



E c 



mi,m2,...,mjv I'rm 



i(*X 



■"ni2 



J2 



mi,m2,...,mjv 



m N 



? mi l r lJ 
V( r 2) 



>m 2 irij 

^m 2 ( r 2) 



y m 2 



JN 



W r iJ 

W( r 2) 



The components of the Slater determinant 
ular orbitals which are usually given as an expansion in "atomic orbitals" 



m .(rA are one-electron molec- 
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Xn- 



(2.17) 



r stays for the Cartesian coordinates (x, y, z) and s is the spin variable 



The collection of coefficients £)... and C... fully characterizes the solution of 
the electronic Schrodinger equation for atoms and molecules. 

The main subject of this course is a discussion of the approximation 
methods for the solution of the Schrodinger equation for the electrons (given 
by the coefficients C... and D,_), which provides the potential for the nuclei 
dynamics V(R). 

Before starting, we have however to translate this problem into a formulation 
suited for computation. Using appropriate basis function it is possible to 
translate 2.10 into a simple linear algebra problem, which can be solved 
using efficient computer software (see for instance the parallel package for 
the solution of linear algebra problem LAPACK: Linear Algebra PACKage 
(www.netlib.org/lapack) 2 ). 



2 LAPACK is written in Fortran77 and provides routines for solving systems of si- 
multaneous linear equations, least-squares solutions of linear systems of equations, eigen- 
value problems, and singular value problems. The associated matrix factorizations (LU, 
Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are related compu- 
tations such as reordering of the Schur factorizations and estimating condition numbers. 
Dense and banded matrices are handled, but not general sparse matrices. In all areas, 
similar functionality is provided for real and complex matrices, in both single and double 
precision. 

The original goal of the LAPACK project was to make the widely used EISPACK and 
LINPACK libraries run efficiently on shared-memory vector and parallel processors. On 
these machines, LINPACK and EISPACK are inefficient because their memory access pat- 
terns disregard the multi-layered memory hierarchies of the machines, thereby spending 
too much time moving data instead of doing useful floating-point operations. LAPACK 
addresses this problem by reorganizing the algorithms to use block matrix operations, 
such as matrix multiplication, in the innermost loops. These block operations can be 
optimized for each architecture to account for the memory hierarchy, and so provide a 
transportable way to achieve high efficiency on diverse modern machines. We use the 
term "transportable" instead of "portable" because, for fastest possible performance, LA- 
PACK requires that highly optimized block matrix operations be already implemented on 
each machine. 

LAPACK routines are written so that as much as possible of the computation is per- 
formed by calls to the Basic Linear Algebra Subprograms (BLAS). While LINPACK and 
EISPACK are based on the vector operation kernels of the Level 1 BLAS, LAPACK was 
designed at the outset to exploit the Level 3 BLAS - a set of specifications for Fortran 
subprograms that do various types of matrix multiplication and the solution of triangular 
systems with multiple right-hand sides. Because of the coarse granularity of the Level 3 
BLAS operations, their use promotes high efficiency on many high-performance comput- 
ers, particularly if specially coded implementations are provided by the manufacturer. 
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2.3 Basis sets, linear algebra and the secular 
equation 



2.3.1 Basis Kets and matrix representation 

Given an Hermitian operator A, its eigenkets (eigenfunctions),|<yC a ) form a 
complete orthonormal set. An arbitrary ket, |a) can be expanded in terms 
of the eigenkets of A. 

\a) = J2 c a\<Pa) (2.18) 

a 

Multiplying with (tp a '\ on the left and using the orthogonality property 
(yv|y?a), we can immediately find the coefficient, 

c a = (<p a \a) (2.19) 

In other words, we have 

|a) = X)l^')W|a), (2.20) 

a' 

which is analogous to an expansion of a vector v in the (real) Euclidean space: 

v = s ^e i (e i -v), (2.21) 

i 

where {ej} form an orthogonal set of unit vectors. 

An important operator is the projection operator A , which acting on a 
ket state \a) gives the component of the ket parallel to \(p a ), 

(MM\) | a) = MM\ a ) = CaM . (2.22) 
Since the sum of all projections of a ket \a) gives back the same ket, 

Y,M(<Pa\ = 1 (2.23) 

where 1 is the identity operator. This representation of the unity operator is 
called the completness relation. 



Highly efficient machine-specific implementations of the BLAS are available for many mod- 
ern high-performance computers. For details of known vendor- or ISV-provided BLAS, 
consult the BLAS FAQ. Alternatively, the user can download ATLAS to automatically 
generate an optimized BLAS library for the architecture. A Fortran77 reference imple- 
mentation of the BLAS is available from netlib; however, its use is discouraged as it will 
not perform as well as a specially tuned implementations. 
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Having specified the base ket, we now show how to represent an operator 
X, by a square matrix. Using the completness relation twice, we can write 
the operator X as 



X = ^2Y1 \<Pa')(<Pa>\X\(p a )(<p a \ 



(2.24) 



There are alltogether N 2 numbers of the form (tp' a \X\(p a ), where N is the 
dimensionality of the ket space. We may arrange them into a N x N square 
matrix where the column and row indices appear as 



'(<^i|X|y?i) ((p!\X\(p2) 

x = I {ip 2 \x\ip 2 ) 



(2.25) 



where the symbol = stands for "is represented by". 

Knowing all (infinite many) matrix elements ((f' a \X\tp a ) of the operator X 
is equivalent to the knowledge of the operator itself (in the same way as 
knowing the 3 componets of a vector in the Euclidean space is sufficient to 
determine its orientation and length). 

In the same way we describe operators by matrices, kets can be described 
by colum vectors, 

/ (<pi\<*)\ 

(<P2\oc) 
(<P3\<X) 



| a) = 



V 



(2.26) 



/ 



and bras as row vectors, 



= {{PM (PM (PM 
= {(<Pi\py faW toW 



(2.27) 



Therefore, the action of an operator on a ket can be represented as a matrix 
multiplication with a vector (link to your linear algebra course). 



2.3.2 Basis functions in quantum chemistry 

In one of the most frequent approximations used in quantum chemistry, the 
complex one-electron or even many-electron molecular wavefunctions are de- 
scribed in basis of atom centered functions. These simplified atomic orbitals 
are often always taken to have the form of sums of Gaussians centered on 
the atoms times a polynomial, Pi, in the electron coordinates relative to that 
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clIjOIXI " 

Xn (r) = ^iVrexp(-ar(|r-R?| 2 ))PKr-R?) ■ (2.28) 
i 

There are conventional sets of these atomic orbitals that are used, that cover 
the polynomials up to a certain order with certain choices of "a"; these are 
called "basis sets" and are given names like "6-31G*" , "TZ2P", and "cc- 
pVQZ". The general procedure is to pick one of these basis sets, and then 
to vary the C's and the D's in 

* el (ri,s 1 ,r 2 , ...,rjv, s N ) = 

= ^2 C , mi,m 2 ,...,m JV |0m 1 (ri, si)0 m2 (r 2 , s 2 ) • • ■ <fim N (r N , s N )\ (2.29) 

mi,m2,...,mjv 

with 

m (r, s) = ^ D mnXn(r) <8> 8 (2.30) 

n 

to try to find an approximate ^ e i ec that solves the Schrodinger equation as 
closely as possible. If your basis set has a very good overlap with the true 
wavefunction, you will be able to achieve good accuracy only varying a few 
C's and D's. 3 



2.3.3 The variational principle and the secular equation 

The variational problem consists in varying the C's and D's to minimize 



E[V elec } = E(C..., £>...) = V e ; e T c| f! c| f c/ (2.31 



elec 


H T elec 


^ elec) 


\^ elec 


^ elect 



For any trial wavefunction,^!/*^, the following inequality holds 



/^f trial 
\ elec 


Helect 


trial \ 
elec 1 


/\Tj trial 
\ elec 


\Tjtrial\ 

elec 1 



El^eiec] < / w nff m 60 • (2-32) 

L eiccj — Uxjtrial \\Titnal\ \ / 



This is called the variational principle. The evaluation of the integral requires 
0(N^ gsis ) operations. (Gaussian functions are used because they allow the 
integrals to be computed analytically.) Typically a basis set might include 
15 atomic orbitals for each atom (except H atoms which do not need so 
many) and you would vary the (15 * A^oms) coefficients D mn . The number 
of possible coefficients C is much larger, something like iVt, as i S raised to the 
-^electrons power, so it is almost always impossible to do anything with a 
complete expansion. Often people don't bother to vary the C's, or only allow 
a small fraction of the C's to vary independently, to reduce the number of 



3 More about the specific basis functions used in computational quantum chemistry will 
follow in Chapter 3. 



32 



parameters. By allowing the C's to vary, you are allowing to account for the 
fact that the different electrons are correlated with each other: when one is 
close to the nucleus the others are likely to be far away. 



2.3.4 Linear variational calculus 

In variational calculus, stationary states of the energy functional are found 
within a subspace of the Hilbert space. An important example is linear vari- 
ational calculus, in which the subspace is spanned by a set of basis vectors 
|S m ),m = 1, ...,M , that we take to be orthonormal. Here we consider the 
case of fixed atomic orbital expansion coefficients (D ) and to-be-optimized 
Slater expansion coefficients (C...) (for example a set of M Slater determi- 
nants, |H m ) = |0 mi (ri)0 m2 (r 2 ) . . .4> mN (r N )\). For a state 

^ e z( r l> r 2, • • • jTat) = ^ Cmi,m2,...,mjvl0mi( r l)0m2( r 2) • • • </W ( r iv) | 

mi,77l2vj m iV 

M 

= ^ y c m |^, m ) (2.33) 

m=l 



the energy functional is given by 



E 



E 



M 



c*c H 



E 



M 
P>9= 



d c p c q °p,q 



(2.34) 



with 

Hp^q — (zip\H e i ec \ziq) (2.35) 

The stationary states follow from the condition that the derivative of this 
functional with respect to the c p vanishes, which leads to 

M 

Y,(H P , q - E 5 M ) c q = 0, for p=l,...,M. (2.36) 

9=1 

Equation (2.36) is an eigenvalue problem which can be written in the matrix 
notation 

HC=EC (2.37) 

This is the Schrodinger equation formulated for a finite, orthonormal basis. 
Although in principle it is possible to use nonlinear parameterizations of the 
wave function, linear parameterizations are used in the large majority of cases 
because of the simplicity of the resulting method, allowing for numerical 
matrix diagonalization techniques. The lowest eigenvalue of (2.37) is 
always higher than or equal to the exact ground state energy, as the 
ground state is the minimal value assumed by the energy functional 
in the full Hilbert space. If we restrict ourselves to a part of this space, 
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then the minimum value of the energy functional must always be higher than 
or equal to the ground state of the full Hilbert space. Including more basis 
functions into our set, the subspace becomes larger, and consequently the 
minimum of the energy functional will decrease (or stay the same). For the 
specific case of linear variational calculus, this result can be generalized to 
higher stationary states: they are always higher than the equivalent solution 
to the full problem, but approximate the latter better with increasing basis 
set size. 

Because the computer time needed for matrix diagonalization scales with the 
third power of the linear matrix size (it is called a 0(M 3 ) process), the basis 
should be kept as small as possible. Therefore, it must be chosen carefully: 
it should be possible to approximate the solutions to the full problem with a 
small number of basis functions 4 . 

In the case in which the basis consists of nonorthonormal basis functions, 

as is often the case in practical calculations, we must reformulate (2.37), 
taking care of the fact that the overlap matrix S, whose elements S pq are 
given by 

Sp, q = (S p |Hq) (2.38) 

is not the unit matrix. This means that in Eq. (2.34) the matrix elements 
S pq of the unit matrix, occurring in the denominator, have to be replaced by 
Spq, and we obtain (for the derivation see the next section) 

HC = £SC. (2.39) 

This looks like an ordinary eigenvalue equation, the only difference being the 
matrix S in the right hand side. It is called a generalized eigenvalue equation 
and there exist computer programs for solving such a problem. 



2.4 Overview of possible approximate solutions 
of the electronic Schrodinger equation 

The most commonly used approximate methods for the solution of the elec- 
tronic molecular Schrodinger equation are: 

- Semi-empirical (MNDO, AMI, PM3, etc.): use a single Slater deter- 
minant (only one C is equal 1 while all the others are set to 0). Vary the 
coefficients D, but just use empirical estimates rather than the true integrals. 
Very cheap, but only accurate for molecule similar to those used to develop 



4 The fact that the basis in (continuous) variational calculus can be chosen so much 
smaller than the number of grid points in a finite difference approach implies that even 
though the latter can be solved using special O(N) methods for sparse systems, they 
are still far less efficient than variational methods with continuous basis functions in most 
cases. This is the reason why, in most electronic structure calculations, variational calculus 
with continuous basis functions is used to solve the Schrodinger equation. 
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the empirical estimates. 

- DFT (B3LYP, BLYP, PW91, etc.): slightly empirical, but much more 
reliable than semi-empirical methods. CPU: cheap, same as HF 0(N 3 ). Er- 
rors ~ 4 kcal/mole (comparable accuracy to MP2 but much cheaper). Pre- 
ferred method for geometries, second derivatives, transition-metal containing 
systems. 

- HF (Hartree-Fock, SCF): only one many-electrons Slater determinant 
is used. Vary the D's. All terms calculated 'ab-initio' within the mean field 
approximation, no empirical parameters. 

CPU: cheap 0(N 3 ) errors ~ 15 kcal/mol. 

- MP2, MP4 (Moller-Plesset, MBPT): Vary the D's first, then set the 
Cs to the values given by perturbation theory (you don't freely vary these 
C's, saving CPU). 

MP2: medium CPU: 0(N 5 ), errors ~ 5 kcal/mol. 

- CI, CISD, QCISD (Configuration Interaction): Vary the coefficients 
D first, freeze them, then vary a lot of the coefficients C. 

Expensive. Not used much anymore, CCSD is preferred. 

- MCSCF, CASSCF: vary a finite set of C's and all the D's simulta- 
neously. Expensive. Good for understanding cases where several electronic 
states have comparable energies. User expertise required to select which C's 
to vary. 

- CAS-PT2: Determine the D's and some C's by CASSCF, then deter- 
mine more C's by perturbation theory. 

Not much more expensive than CASSCF. Sometimes very good, but not reli- 
able. 

- MRCI (multi reference CI): Determine the D's and some C's by CASSCF 
or MCSCF, freeze these, then allow many of the C's to vary. 

Super expensive. Very high accuracy for small systems. 

- CCSD, CCSD(T) (Coupled Cluster): Vary the D's, fix them, then 
vary a lot of the C's, but constraining certain relationships between the C's. 
This allows you to effectively use a longer expansion without increasing the 
number of adjustable parameters so much. The constraints force the solution 
to be "size-consistent", i.e. two molecules calculated simultaneously have ex- 
actly the same energy as two molecules calculated separately. 
Expensive. Often very accurate. 

- Extrapolations ("Composite Methods"): G2, G3, CBS-q, CBS-Q, 
CBS-QB3, CBS-RAD Run a series of the above calculations with different 
size basis sets, following some recipe. The results from all these calculations 
are extrapolated to an estimate of the true potential V(R). These methods 
give excellent accuracy in less CPU time than CCSD or MRCI. However, 
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the multiple steps involved provide many opportunities for something to go 
wrong. 

Accuracy: usually 1-2 kcal/mol. 

Some Practical Warnings 

1) The optimization (SCF/HF/DFT/CASSCF/MRSCF) problem required 
to solve for the D's is nonlinear and has multiple solutions, only one of 
which is the one you want (usually you want the lowest energy solution). 
So you may end up converging to a wavefunction which is qualitatively 
incorrect, perhaps it corresponds to an electronically excited state. 

2) Most of the quantum chemistry methods have problems (convergence, 
accuracy) with systems where there are low-lying electronic states (close 
to the ground state). In these cases, sometimes the numbers computed 
are completely nuts, other times they are subtly wrong. This is par- 
ticularly a problem for transition states and where there are several 
lone pair electrons in the system. If you must study these systems, get 
expert assistance. 

3) Many molecules have multiple geometrical conformations (local minima 
in V(R)), and sometimes there are multiple saddle points that might 
be confused with the transition state (TS). Look at your structures, 
if they are not what you expected, investigate. Also, it is worth some 
effort to make sure your initial guess at the molecular geometry is quite 
good, otherwise the geometry-optimization algorithm may get lost and 
waste a lot of CPU time to no avail. If you are having troubles, you 
can constrain some of the coordinates to make things easier for the 
optimizer. 

4) For radicals and other open-shell systems, compare your computed so- 
lutions (S 2 ) with the theoretical value S(S + 1). If your number is 
way off, chances are you have other problems as well. Sometimes you 
can use "restricted" methods like ROHF and RMP2, or spin-projection 
methods to fix this "spin-contamination" problem. 

5) Every method runs into problems sometimes, and sometimes they are 
quite subtle. It is a good idea to double check your calculation with 
another calculation done using a very different method. If they both 
agree you can be pretty confident that your result is real. 
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Chapter 3 



Basis functions in quantum 
chemistry 



This chapter is adapted from Chapter 5 of Jensen's book: F. Jensen, 'Intro- 
duction to Computational Chemistry', Wiley. 

In the derivation in the previous chapter, we have introduced the concept 
of basis function for the expansion of the one-electron molecular orbitals used 
for the generation of the many-electrons wave functions (Slater determinants 
or linear combination of Slater determinants). 
There we derived the following expansion (eq. 2.17): 



(where \n is an atom centered basis function and the spin dependent part of 
the wavefunctions is left out). 

In this chapter, we introduce the different basis functions, \ n commonly used 
in computational quantum chemistry. 



Finiteness of Basis Sets: Approximations 

One of the approximations inherent in essentially all ab initio methods is 
the introduction of a finite basis set. Expanding an unknown function, such 
as a molecular orbital, in a set of known functions is not an approximation, 
if the basis is complete. However, a complete basis means that an infinite 
number of functions must be used, which is impossible in actual calculations. 
An unknown MO can be thought of as a function in the infinite coordinate 
system spanned by the complete basis set. When a finite basis is used, only 
the components of the MO along those coordinate axes corresponding to 
the selected basis can be represented. The smaller the basis, the poorer the 
representation. The type of basis functions used also influences the accuracy. 




(3.1) 



n 
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The better a single basis function is able to reproduce the unknown function, 
the fewer basis functions necessary for achieving a given level of accuracy. 
Knowing that the computational effort of ab initio methods scales formally 
at least as M 4 , it is of course of prime importance to make the basis set as 
small as possible without compromising the accuracy. 



3.1 Slater and Gaussian Type Orbitals 

There are two types of basis functions (also called Atomic Orbitals, AO, al- 
though in general they are not solutions to an atomic Schrodinger equation) 
commonly used in electronic structure calculations: Slater Type Orbitals 
(STO) and Gaussian Type Orbitals (GTO). 

A procedure that has come into wide use is to fit a Slater-type orbital (STO) 
to a linear combination of n — 1, 2, 3, . . . primitive Gaussian functions. This 
is the STO-nG procedure. In particular, STO-3G basis sets are often used 
in polyatomic calculations, in preference to evaluating integrals with Slater 
functions. 




Radius (a u.) Radius (o.u.) 

Figure 3.1: Comperison of Slater function with Gaussian function: 
a) least squares fit of a Is Slater function (£ = 1.0) by a single STO- 
1G Is Gaussian function; b) comparison of the corresponding radial 
distribution functions (47rr 2 |</(>i s (r)| 2 ). 



1. Slater type orbitals have the functional form 

Xc,n,i,m(r, 6 t tp)=N Y ltm (6, tp) r"" 1 e~< T (3.2) 

N is a normalization constant and Y^ m are the usual spherical harmonic 
functions. The exponential dependence on the distance between the nucleus 
and the electron mirrors the exact decay behavior of the orbitals for the 
hydrogen atom. However, since STOs do not have any radial nodes, nodes in 
the radial part are introduced by making linear combinations of STOs. The 
exponential dependence ensures a fairly rapid convergence with increasing 
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number of functions, however, the calculation of three- and four-centre two 
electron integrals cannot be performed analytically. STOs are primarily used 
for atomic and diatomic systems where high accuracy is required, and in semi- 
empirical methods where all three- and four- center integrals are neglected. 




■0 .5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 
Radius (a.u.) 



Figure 3.2: Comparison of the quality of the least-squares fits of 
a Is Slater function (C = 1.0) obtained at the STO-1G, STO-2G, 
and STO-3G levels. 



2. Gaussian type orbitals can be written in terms of polar or cartesian 
coordinates 

X(,n,i,m(r, 0,<P) = N YiM <P) r2n ~ 2 ~ l e " ?r2 ( 3 - 3 ) 
X(,i x ,i v ,h 0, V,z) = N x x * y l y y l * (3.4) 

where the sum of l x ,l y and l z determines the type of orbital (for example 
l x + ly + lz — 1 is a p-orbital) 1 . 



1 Although a GTO appears similar in the two sets of coordinates, there is a subtle dif- 
ference. A d-type GTO written in terms of the spherical functions has five components 
(^2,2,^2,1,^2,0)^2,-1,^2,-2)7 but there appear to be six components in the Cartesian co- 
ordinates (x2, ?/2, ^2, xy, xz, yz). The latter six functions, however, may be transformed to 
the five spherical d- functions and one additional s- function (x 2 +y 2 + z 2 ). Similarly, there 
are 10 Cartesian "f-functions" which may be transformed into seven spherical f- functions 
and one set of spherical p-functions. Modern programs for evaluating two-electron inte- 
grals are geared to Cartesian coordinates, and they generate pure spherical d-functions 
by transforming the six Cartesian components to the five spherical functions. When only 
one d-function is present per atom the saving by removing the extra s-function is small, 
but if many d-functions and/or higher angular moment functions (f-, g-, h- etc. functions) 
are present, the saving can be substantial. Furthermore, the use of only the spherical 
components reduces the problems of linear dependence for large basis sets, as discussed 
below. 
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3. Comparison between STO and GTO 

i. The r 2 dependence in the exponent makes the GTOs inferior to the 
STOs in two aspects. At the nucleus the GTO has zero slope, in con- 
trast to the STO which has a "cusp" (discontinuous derivative), and 
GTOs have problems representing the proper behavior near the nu- 
cleus. 

ii. The other problem is that the GTO falls off too rapidly far from the 
nucleus compared with an STO, and the "tail" of the wave function is 
consequently represented poorly. 

iii. Both STOs and GTOs can be chosen to form a complete basis, but the 
above considerations indicate that more GTOs are necessary for achiev- 
ing a certain accuracy compared with STOs. A rough guideline says 
that three times as many GTOs as STOs are required for reaching a 
given level of accuracy. The increase in number of basis functions, how- 
ever, is more than compensated for by the ease by which the required 
integrals can be calculated. In terms of computational efficiency, GTOs 
are therefore preferred, and used almost universally as basis functions 
in electronic structure calculations. 



3.2 Classification of Basis Sets 

Having decided on the type of basis function (STO/GTO) and their location 
(nuclei), the most important factor is the number of functions to be used. 
The smallest number of functions possible is a minimum basis set. Only 
enough atomic orbital functions are employed to contain all the electrons of 
the neutral atom(s). 

3.2.1 Minimum basis sets. Examples 

For hydrogen (and helium) this means a single s-function. For the first row 
in the periodic table it means two s-functions (Is and 2s) and one set of 
p-functions (2p x , 2p y and 2p z ). Lithium and beryllium formally only require 
two s-functions, but a set of p-functions is usually also added. For the second 
row elements, three s-functions (Is, 2s and 3s) and two sets of p-functions 
(2p and 3p) are used. 



3.2.2 Improvements 

1. The first improvement in the basis sets is a doubling of all basis func- 
tions, producing a Double Zeta (DZ) type basis. The term zeta stems 
from the fact that the exponent of STO basis functions is often denoted 
by the greek letter (. 

A DZ basis thus employs two s-functions for hydrogen (Is and Is'), four 
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s-functions (Is, Is', 2s and 2s') and two p-functions (2p and 2p') for first 
row elements, and six s-functions and four p-functions for second row 
elements. Doubling the number of basis functions allows for a much 
better description of the fact that the electron distribution in molecules 
can differ significantly from the one in the atoms and the chemical bond 
may introduce directionalities which can not be described by a minimal 
basis. 

The chemical bonding occurs between valence orbitals. Doubling the 
ls-functions in for example carbon allows for a better description of 
the ls-electrons. However, the Is orbital is essentially independent of 
the chemical environment, being very close to the atomic case. A vari- 
ation of the DZ type basis only doubles the number of valence orbitals, 
producing a split valence basis. 2 . 



2. The next step up in basis set size is a Triple Zeta (TZ) basis. Such a 
basis contains three times as many functions as the minimum basis, i.e. 
six s-functions and three p-functions for the first row elements. Some 
of the core orbitals may again be saved by only splitting the valence, 
producing a triple zeta split valence basis set. The names Quadruple 
Zeta (QZ) and Quintuple Zeta (5Z, not QZ) for the next levels of 
basis sets are also used, but large sets are often given explicitly in terms 
of the number of basis functions of each type. 



3. In most cases higher angular momentum functions are also important, 
these are denoted polarization functions. Consider for example a C- 
H bond which is primarily described by the hydrogen s-orbital(s) and 
the carbon s- and p 2 -orbitals. It is clear that the electron distribution 
along the bond will be different than that perpendicular to the bond. If 
only s-functions are present on the hydrogen, this cannot be described. 
However, if a set of p-orbitals is added to the hydrogen, the p com- 
ponent can be used for improving the description of the H-C bond. 
The p-orbital introduces a polarization of the s-orbital(s). Similarly, 
d-orbitals can be used for polarizing p-orbitals, f-orbitals for polarizing 
d-orbitals etc. Once a p-orbital has been added to a hydrogen s-orbital, 
it may be argued that the p-orbital now should be polarized by adding a 
d-orbital, which should be polarized by an f-orbital, etc. For single de- 
terminant wave functions, where electron correlation is not considered, 
the first set of polarization functions (i.e. p-functions for hydrogen and 
d-functions for heavy atoms) is by far the most important, and will in 
general describe all the important charge polarization effects. 
Adding a single set of polarization functions (p-functions on hydrogens 
and d-functions on heavy atoms) to the DZ basis forms a Double Zeta 



2 In actual calculations a doubling of the core orbitals would rarely be considered, and 
the term DZ basis is also used for split valence basis sets (or sometimes VDZ, for valence 
double zeta) 



41 



plus Polarization (DZP) type basis 3 . Similarly to the sp-basis sets, 
multiple sets of polarization functions with different exponents may be 
added. If two sets of polarization functions are added to a TZ sp-basis, 
a Triple Zeta plus Double Polarization (TZ2P) type basis is obtained. 
For larger basis sets with many polarization functions the explicit com- 
position in terms of number and types of functions is usually given. 
At the HF level there is usually little gained by expanding the basis 
set beyond TZ2P, and even a DZP type basis set usually gives "good" 
results (compared to the HF limit). 



3.3 Basis set balance 

In principle many sets of polarization functions may be added to a small 
sp-basis. This is not a good idea. If an insufficient number of sp-functions 
bas been chosen for describing the fundamental electron distribution, the 
optimization procedure used in obtaining the wave function (and possibly 
also the geometry) may try to compensate for inadequacies in the sp-basis 
by using higher angular momentum functions, producing artefacts. A rule of 
thumb says that the number of functions of a given type should at most be 
one less than the type with one lower angular momentum. A 3s2pld basis is 
balanced, but a 3s2p2d2flg basis is too heavily polarized. 
Another aspect of basis set balance is the occasional use of mixed basis 
sets, for example a DZP quality on the atoms in the "interesting" part of the 
molecule and a minimum basis for the "spectator" atoms. Another example 
would be addition of polarization functions for only a few hydrogens which 
are located "near" the reactive part of the system. For a large molecule this 
may lead to a substantial saving in the number of basis functions. It should 
be noted that this may bias the results and can create artefacts. For example, 
a calculation on the H 2 molecule with a minimum basis at one end and a DZ 
basis at the other end will predict that H2 has a dipole moment, since the 
variational principle will preferentially place the electrons near the center 
with the most basis functions. The majority of calculations are therefore 
performed with basis sets of the same quality (minimum, DZP, TZ2P, . . .) 
on all atoms, possibly cutting polarization and/or diffuse (small exponent) 
functions on hydrogens. 

Except for very small systems it is impractical to saturate the basis set so 
that the absolute error in the energy is reduced below chemical accuracy, for 



3 There is a variation where polarization functions are only added to non- hydrogen 
atoms. This does not mean that polarization functions are not important on hydrogens. 
However, hydrogens often have a "passive" role, sitting at the end of bonds which does 
not take an active part in the property of interest. The errors introduced by not including 
hydrogen polarization functions are often rather constant and, as the interest is usually in 
energy differences, they tend to cancel out. As hydrogens often account for a large number 
of atoms in the system, a saving of three basis functions for each hydrogen is significant. 
If hydrogens play an important role in the property of interest, it is of course not a good 
idea to neglect polarization functions on hydrogens. 
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example 1 kcal/ mol. The important point in choosing a balanced basis set is 
to keep the error as constant as possible. The use of mixed basis sets should 
therefore only be done after careful consideration. Furthermore, the use of 
small basis sets for systems containing elements with substantially different 
numbers of valence electrons (like LiF) may produce artefacts. 

3.4 How do we choose the exponents in the ba- 
sis functions? 

The values for s- and p-functions are typically determined by performing 
variational HF calculations for atoms, using the exponents as variational pa- 
rameters. The exponent values which give the lowest energy are the "best", 
at least for the atom. In some cases the optimum exponents are chosen on 
the basis of minimizing the energy of a wave function which includes electron 
correlation. The HF procedure cannot be used for determining exponents of 
polarization functions for atoms. By definition these functions are unoccu- 
pied in atoms, and therefore make no contribution to the energy. Suitable 
polarization exponents may be chosen by performing variational calculations 
on molecular systems (where the HF energy does depend on polarization 
functions) or on atoms with correlated wave functions. Since the main func- 
tion of higher angular momentum functions is to recover electron correlation, 
the latter approach is usually preferred. Often only the optimum exponent 
is determined for a single polarization function, and multiple polarization 
functions are generated by splitting the exponents symmetrically around the 
optimum value for a single function. The splitting factor is typically taken in 
the range 2-4. For example if a single d-function for carbon has an exponent 
value of 0.8, two polarization functions may be assigned with exponents of 
0.4 and 1.6 (splitting factor of 4). 



3.5 Contracted Basis functions 

One disadvantage of all energy optimized basis sets is the fact that they 
primarily depend on the wave function in the region of the inner shell elec- 
trons. The ls-electrons account for a large part of the total energy, and 
minimizing the energy will tend to make the basis set optimal for the core 
electrons, and less than optimal for the valence electrons. However, chemistry 
is mainly dependent on the valence electrons. Furthermore, many proper- 
ties (for example polarizability) depend mainly on the wave function "tail" 
(far from the nucleus), which energetically is unimportant. An energy opti- 
mized basis set which gives a good description of the outer part of the wave 
function needs therefore to be very large, with the majority of the functions 
being used to describe the ls-electrons with an accuracy comparable to that 
for the outer electrons in an energetic sense. This is not the most efficient 
way of designing basis sets for describing the outer part of the wave function. 
Instead energy optimized basis sets are usually augmented explicitly with 
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diffuse functions (basis functions with small exponents). Diffuse functions 
are needed whenever loosely bound electrons are present (for example in an- 
ions or excited states) or when the property of interest is dependent on the 
wave function tail (for example polarizability). 

The fact that many basis functions go into describing the energetically im- 
portant, but chemically unimportant, core electrons is the foundation for 
contracted basis sets. 

An example. The carbon atom Consider for example a basis set consist- 
ing of 10 s-mnctions (and some p-functions) for carbon. Having optimized these 
10 exponents by variational calculations on the carbon atom, maybe six of the 10 
functions are found primarily to be used for describing the Is orbital, and two of the 
four remaining describe the "inner" part of the 2s-orbital. The important chemical 
region is the outer valence. Out of the 10 functions, only two are actually used 
for describing the chemically interesting phenomena. Considering that the com- 
putational cost increases as the fourth power (or higher) of the number of basis 
functions, this is very inefficient. As the core orbitals change very little depending 
on the chemical bonding situation, the MO expansion coefficients in front of these 
inner basis functions also change very little. The majority of the computational 
effort is therefore spent describing the chemically uninteresting part of the wave 
function, which furthermore is almost constant. Consider now making the varia- 
tional coefficients in front of the inner basis functions constant, i.e. they are no 
longer parameters to be determined by the variational principle. The ls-orbital is 
thus described by a fixed linear combination of say six basis functions. Similarly 
the remaining four basis functions may be contracted into only two functions, for 
example by fixing the coefficient in front of the inner three functions. In doing this 
the number of basis functions to be handled by the variational procedure has been 
reduced from 10 to three. 

Combining the full set of basis functions, known as the primitive GTOs 
(PGTOs), into a smaller set of functions by forming fixed linear combina- 
tions is known as basis set contraction, and the resulting functions are 
called contracted GTOs (CGTOs) 

k 

X (CGTO) = ^a lXi (PGTO) (3.5) 

i 

The previously introduced acronyms DZP, TZ2P etc., refer to the number 
of contracted basis functions. Contraction is especially useful for orbitals 
describing the inner (core) electrons, since they require a relatively large 
number of functions for representing the wave function cusp near the nucleus, 
and furthermore are largely independent of the environment. Contracting 
a basis set will always increase the energy, since it is a restriction of the 
number of variational parameters, and makes the basis set less flexible, but 
will also reduce the computational cost significantly. The decision is thus how 
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much loss in accuracy is acceptable compared to the gain in computational 
efficiency. 

3.5.1 The degree of contraction 

The degree of contraction is the number of PGTOs entering the CGTO, 
typically varying between 1 and 10. The specification of a basis set in terms 
of primitive and contracted functions is given by the notation 

(10s4pld/4slp) — > [3s2pld/2slp] . (3.6) 

The basis in parentheses is the number of primitives with heavy atoms (first 
row elements) before the slash and hydrogen after. The basis in the square 
brackets is the number of contracted functions. Note that this does not tell 
how the contraction is done, it only indicates the size of the final basis (and 
thereby the size of the variational problem in HF calculations). 

3.6 Example of Contracted Basis Sets; Pople 
Style Basis Sets 

There are many different contracted basis sets available in the literature or 
built into programs, and the average user usually only needs to select a suit- 
able quality basis for the calculation. For short description of some basis 
sets which often are used in routine calculations (see for instance the book of 
Frank Jensen, Introduction to Computational Chemistry, Wiley, 2002. Chap- 
ter 5). 



STO-nG basis sets n PGTOs fitted to a 1 STO. This is a minimum 
type basis where the exponents of the PGTO are determined by fitting to 
the STO, rather than optimizing them by a variational procedure. Although 
basis sets with n = 2 — 6 have been derived, it has been found that using 
more than three PGTOs to represent the STO gives little improvement, and 
the STO-3G basis is a widely used minimum basis. This type of basis set has 
been determined for many elements of the periodic table. The designation 
of the carbon/hydrogen STO-3G basis is (6s3p/3s) — ► [2slp/ls]. 

k-nlmG basis sets These basis sets have been designed by Pople and co- 
workers, and are of the split valence type, with the k in front of the dash 
indicating how many PGTOs are used for representing the core orbitals. The 
nlm after the dash indicate both how many functions the valence orbitals 
are split into, and how many PGTOs are used for their representation. Two 
values (e.g. nl) indicate a split valence, while three values (e.g. nlm) indicate 
a triple split valence. The values before the G (for Gaussian) indicate the s- 
and p-functions in the basis; the polarization functions are placed after the 
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G. This type of basis sets has the further restriction that the same exponent 
is used for both the s- and p-functions in the valence. This increases the com- 
putational efficiency, but of course decreases the flexibility of the basis set. 
The exponents in the PGTO have been optimized by variational procedures. 

3-21G This is a split valence basis, where the core orbitals are a contraction 
of three PGTOs, the inner part of the valence orbitals is a contraction of two 
PGTOs and the outer part of the valence is represented by one PGTO. The 
designation of the carbon/hydrogen 3-21G basis is (6s3p/3s) — > [3s2p/2s]. 
Note that the 3-2 1G basis contains the same number of primitive GTOs as 
the STO-3G, however, it is much more flexible as there are twice as many 
valence functions which can combine freely to make MOs. 

6-31G This is also a split valence basis, where the core orbitals are a con- 
traction of six PGTOs, the inner part of the valence orbitals is a contraction 
of three PGTOs and the outer part of the valence represented by one PGTO. 
The designation of the carbon/hydrogen 6-31G basis is (10s4p/4s) — > [3s2p/2s]. 
In terms of contracted basis functions it contains the same number as 3-2lG, 
but the representation of each functions is better since more PGTOs are 
used. 

6-311G This is a triple zeta split valence basis, where the core orbitals 
are a contraction of six PGTOs and the valence split into three functions, 
represented by three, one, and one PGTOs, respectively. 

To each of these basis sets one can add diffuse and/or polarization 
functions. 

• Diffuse functions are normally s- and p-functions and consequently go 
before the G. They are denoted by + or ++, with the first + indicating 
one set of diffuse s- and p-functions on heavy atoms, and the second + 
indicating that a diffuse s-function is also added to hydrogens. The ar- 
guments for adding only diffuse functions on non-hydrogen atoms is the 
same as that for adding only polarization functions on non-hydrogens. 

• Polarization functions are indicated after the G, with a separate 
designation for heavy atoms and hydrogens. The 6-31+G(d) is a split 
valence basis with one set of diffuse sp-functions on heavy atoms only 
and a single d-type polarization function on heavy atoms. 

A 6-311 + + G(2df,2pd) is similarly a triple zeta split valence with ad- 
ditional diffuse sp-functions, and two d- and one f-functions on heavy 
atoms and diffuse s- and two p- and one d-functions on hydrogens. The 
largest standard Pople style basis set is 6-311 ++G(3df, 3pd). These 
types of basis sets have been derived for hydrogen and the first row 
elements, and same of the basis sets have also been derived for second 
and higher row elements. 

If only one set of polarization functions is used, an alternative notation 
in terms of * is also widely used. The 6-31 G* basis is identical to 6- 
31G(d), and 6-31G** is identical to 6-31G(d,p). A special note should 
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be made for the 3-21G* basis. The 3-21G basis is basicly too small to 
support polarization functions (it becomes unbalanced). However, the 
3-21 G basis by itself performs poorly for hypervalent molecules, such as 
sulfoxides and sulfones. This can be substantially improved by adding 
a set of d-functions. The 3-21 G* basis has only d-functions on second 
row elements (it is sometimes denoted 3-21 G(*) to indicate this), and 
should not be considered a polarized basis. Rather, the addition of a 
set of d-functions should be considered an ad hoc repair of a known 
flaw. 
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Chapter 4 



An Introduction to Hartree Fock 
Theory 

Adapted from C. D. Sherrill's notes: "An Introduction to Hartree-Fock Molec- 
ular Orbital Theory". 



4.1 Introduction 



Hartree-Fock theory is fundamental to much of electronic structure theory. 
It is the basis of molecular orbital (MO) theory, which posits that each elec- 
tron's motion can be described by a single-particle function (orbital) which 
does not depend explicitly on the instantaneous motions of the other elec- 
trons. Some of you have probably learned about (and maybe even solved 
problems with) Hiickel MO theory, which takes Hartree-Fock MO theory 
as an implicit foundation and throws away most of the terms to make it 
tractable for simple calculations. The ubiquity of orbital concepts in chem- 
istry is a testimony to the predictive power and intuitive appeal of Hartree- 
Fock MO theory. However, it is important to remember that these orbitals 
are mathematical constructs which only approximate reality. Only for the 
hydrogen atom (or other one-electron systems, like He + ) are orbitals exact 
eigenfunctions of the full electronic Hamiltonian. As long as we are content 
to consider molecules near their equilibrium geometry, Hartree-Fock theory 
often provides a good starting point for more elaborate theoretical meth- 
ods which are better approximations to the electronic Schrodinger equation 
(e.g., many-body perturbation theory, single- reference configuration interac- 
tion). So, how do we calculate molecular orbitals using Hartree-Fock theory? 
This is the subject of these notes; we will explain Hartree-Fock theory at an 
introductory level. 
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4.2 What Problem Are We Solving? 



It is always important to remember the context of a theory. Hartree-Fock the- 
ory was developed to solve the electronic Schrodinger equation that results 
from the time-independent Schrodinger equation after invoking the Born- 
Oppenheimer approximation. In atomic units, and with r denoting elec- 
tronic and R denoting nuclear degrees of freedom, the electronic Schrodinger 
equation is 

t 7,i 1 J 11 7>J 1 J J| i>j 1 1 11 

(4.1) 

or, in a more compact notation, 

[f e (r) + V^(r, R) + VWtv(R) + V ee (r)] tt(r, R) = £ e ^(r, R) (4.2) 

Recall from the Born-Oppenheimer approximation that E e i (plus or minus 

Vjvtv(R), which we include here) will give us the potential energy experi- 
enced by the nuclei. In other words, E e i(R) the potential energy surface 
(from which we can get, for example, the equilibrium geometry vibrational 
frequencies). That's one good reason why we want to solve the electronic 
Schrodinger equation. The other is that the electronic wavefunction ^(r, R) 
contains lots of useful inform about molecular properties such as dipole (and 
multipole) moments, polarizability, etc. 

4.3 The many-electron wavefunction: the Slater 
determinant 

The basic idea of Hartree-Fock theory is as follows. We know how to solve 
the electronic problem for the simplest atom, hydrogen, which has only one 
electron. We imagine that perhaps if we added another electron to hydrogen, 
to obtain H~, then maybe it might be reasonable to start off pretending that 

the electrons don't interact with each other (i.e., that V ee = 0). If that was 
true, then the Hamiltonian would be separable, and the total electronic wave- 
function \&(ri,r2) describing the motions of the two electrons would just be 
the product of two hydrogen atom wavefunctions (orbitals), \I/#(ri) ^h{^2) 
(you should be able to prove this easily). 

However, we have already shown in the previous chapter that a correct 
many electron wavefunction must fulfill the antisymmetry principle together 
with the principle of indistinguishability of the electrons. The solution of this 
problem is the Slater determinant made of one-electron molecular orbitals 
(MOs). In the following we will derive a theory, the Hartree Fock theory, for 
the calculation of such single electron MOs. 
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4.4 Simplified Notation for the Hamiltonian 



Hartree-Fock theory is based on the assumption that the many-electron wave- 
function of the system can be described by a single Slater determinant made 
of one-electron molecular orbitals. Let's re-examine the Hamiltonian to make 
it look as simple as possible. In the process, we will bury some complexity 
that would have to be taken care of later (in the evaluation of integrals) . 

We will define a one-electron operator h as follows 
and a two-electron operator v(i,j) as 



v(i,j) = , _ , (4.4) 

I r i r j I 

Now we can write the electronic Hamiltonian much more simply, as 

H el = ^h{i) + + V NN (4.5) 



i<3 



Since Vnn is just a constant for the fixed set of nuclear coordinates {R}, it 
doesn't change the eigenfunctions, and only shifts the eigenvalues. 



4.5 Energy Expression 

Now that we have a form for the wavefunction and a simplified notation for 
the Hamiltonian, we have a good starting point to tackle the problem. Still, 
how do we obtain the molecular orbitals? 

We state that the Hartree-Fock wavefunction will have the form of a 
Slater determinant, and that the energy will be given by the usual quantum 
mechanical expression (assuming the wavefunction is normalized): 

E el = (V\H el \V). (4.6) 

where \1/ = ^(i - !, r 2 , . . . , rjv) is a Slater determinant, 2 , . . . , 

We can employ the variational theorem, which states that the energy is 
always an upper bound to the true energy. Hence, we can obtain better ap- 
proximate wavefunctions \l/ by varying their parameters until we minimize 
the energy within the given functional space. Hence, the correct molecular 
orbitals in the Slater determinant, {4>i}fL\, are those which minimize the elec- 
tronic energy E e i\ The molecular orbitals can be obtained numerically using 
integration over a grid, or (much more commonly) as a linear combination 
of a set of given basis functions (so-called "atomic orbital" basis functions, 
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usually atom-centered Gaussian type functions). 

Now, using some tricks we don't have time to get into, we can re-write 
the Hartree-Fock energy E e i in terms of integrals of the one- and two-electron 
operators: 

Ehf = J2(i\h\i) + - X>1jj] - (4.7) 

i ij 

where the one electron integral is (xj = (r», Sj)) 1 

(i\h\j) = j dx 1 0*(x 1 )/ i (r 1 )^(x 1 ) (4.8) 
and a two-electron integral (Chemists' notation) is 

[ij\kl] = [ dx 1 dx. 2 <f>*{xi)<f> j (x 1 )— fa k {x 2 )fa(x 2 ). (4.9) 

There exist efficient computer algorithms for computing such one- and two- 
electron integrals. 

4.6 The Hartree-Fock Equations 

Again, the Hartree-Fock method seeks to approximately solve the electronic 
Schrodinger equation, and it assumes that the wavefunction can be approxi- 
mated by a single Slater determinant made up of one spin orbital per electron. 
According to the variational theorem we know that the Slater determinant 
with the lowest energy is as close as we can get to the true wavefunction for 
the assumed functional form of a single Slater determinant. The Hartree- 
Fock method determines the set of spin orbitals which minimize the energy 
and give us this "best single determinant." 

So, we need to minimize the Hartree-Fock energy expression with respect 
to changes in the orbitals fa — ► fa + 5 fa. We have also been assuming that 
the orbitals <ft are orthonormal, and we want to ensure that our variational 
procedure leaves them orthonormal. We can accomplish this by Lagrange's 
method of undetermined multipliers, where we employ a functional £ defined 
as 

C[{fa}] = E HF [{fa}] »,((/./) - <y (4.10) 

ij 

where are the undetermined Lagrange multipliers and is the overlap 
between spin orbitals i and j, i.e., 
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Setting the first variation 8C = 0, and working through some algebra, we 
eventually arrive at the Hartree-Fock equations defining the orbitals: 



h(xx)<f>j(xi) + J dx 2 1 fa (x 2 



-E 



dx 2 0*(x 2 )0 i (x 2 )r 12 1 



[xxj = ei0i(xi, 



(4.12) 
(4.13) 



where 6{ is the energy eigenvalue associated with orbital fa. 

The Hartree-Fock equations can be solved numerically (exact Hartree- 
Fock), or they can be solved in the space spanned by a set of basis functions 
(Hartree-Fock- Roothan equations). In either case, note that the solutions 
depend on the orbitals. Hence, we need to guess some initial orbitals and 
then refine our guesses iteratively. For this reason, Hartree-Fock is called a 
self-consistent-field (SCF) approach. 

The first term above in square brackets, 



E 



^ x 2|0j(x 2 )| 2 r 12 1 



Xl 



(4.14) 



(r 12 = |ri — r 2 |) gives the Coulomb interaction of an electron in spin orbital 
fa with the average charge distribution of the other electrons. Here we see in 
what sense Hartree-Fock is a "mean field" theory. This is called the Coulomb 
term, and it is convenient to define a Coulomb operator as 



tfac 2 |^(x2)| 2 r 12 1 , 



(4.15) 



which gives the average local potential at point Xi due to the charge distri- 
bution from the electron in orbital <pj. 

The other term in brackets in eq.(20) is harder to explain and does not 
have a simple classical analog. It arises from the antisymmetry requirement 
of the wavefunction. It looks much like the Coulomb term, except that it 
switches or exchanges spin orbitals fa and <pj . Hence, it is called the exchange 
term: 



E 



rfx 2 0*(x 2 )0 i (x 2 )r 12 1 



•i x i ■ 



(4.16) 



We can define an exchange operator in terms of its action on an arbitrary 
spin orbital fa 



/C j (xi)0 i (xi 



dx 2 0*(x 2 )r 12 Vi(x 2 ) 



(4.17) 



In terms of these Coulomb and exchange operators, the Hartree-Fock equa- 
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tions become considerably more compact. 



XiJ = 6i0i(Xi 



(4.18) 



Perhaps now it is more clear that the Hartree-Fock equations are eigenvalue 
equations. If we realize that 



[j;(x 1 )-/c l (x 1 )]^(x 1 ) = o, 



(4.19) 



then it becomes clear that we can remove the restrictions j '• ^ i in the sum- 
mations, and we can introduce a new operator, the Fock operator, as 



/(x 1 ) = Mxi) + ^^(x 1 )-/C,(x 1 ). 



(4.20) 



And now the Hartree-Fock equations are just 

/(Xi)^i(Xi) = €i0i(Xi) 



(4.21) 



4.6.1 Matrix representation of the Hartree-Fock equa- 
tion: the Roothaan equation 



Introducing a basis set transforms the Hartree-Fock equations into the Roothaan 
equations. Denoting the atomic orbital basis functions as \i we have the ex- 
pansion 

K 

11=1 

for each spin orbital i. This leads to 

/( Xl ) D *X v {xi) = ti A^xi). ( 4 -23) 

V V 

Left multiplying by x^(xi) and integrating yields a matrix equation 

^D ui / dx 1 x*(x 1 )/(xi)x„(xi) = ei^Dvi / rfxix*(xi)x^(xi). (4.24) 
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This can be simplified by introducing the matrix element notation 



S^u = J rfxix*(xi)x i ,(xi), (4.25) 
F, u = I dx lX :(x 1 )/(x 1 )x„(x 1 ). (4.26) 



Now the Hartree-Fock-Roothaan equations can be written in matrix form as 

F^D vi = €i S^Dvi (4.27) 

V V 

or even more simply as matrices 

FD = SDe (4.28) 

where e is a diagonal matrix of the orbital energies e^. This is like an eigen- 
value equation except for the overlap matrix S. One performs a transforma- 
tion of basis to go to an orthogonal basis to make S vanish. Then it's just 
a matter of solving an eigenvalue equation (or, equivalently, diagonalizing 
F!). Well, not quite. Since F depends on it's own solution (through the 
orbitals), the process must be done iteratively. This is why the solution of 
the Hartree-Fock-Roothaan equations are often called the self- consistent- field 
procedure. 



Computational Aspects: Variational Optimiza- 
tion of Orbitals 

The variational theorem states that for a time-independent Hamiltonian op- 
erator, any trial wavefunction will have an energy expectation value that is 
greater than or equal to the true ground state wavefunction corresponding 
to the given Hamiltonian. 

E = i ^ir - Eo (429) 

Because of this, the Hartree-Fock energy is an upper bound to the true ground 
state energy of a given molecule. In the context of the Hartree-Fock method, 
the best possible solution is at the Hartree-Fock limit, i.e. the limit of the 
Hartree-Fock energy as the basis set approaches completeness. (The other 
is the full-CI limit, where the last two approximations of the Hartree-Fock 
theory as described above are completely undone. It is only when both limits 
are attained that the exact solution is obtained.) 

The starting point for the Hartree-Fock method is a set of approximate 
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one-electron wavef unctions known as orbitals. For an atomic calculation, 
these are typically the orbitals for a hydrogenic atom (an atom with only 
one electron, but the appropriate nuclear charge). For a molecular or crys- 
talline calculation, the initial approximate one-electron wavefunctions are 
typically a linear combination of atomic orbitals (LCAO). 

The orbitals above only account for the presence of other electrons in an 
average manner. In the Hartree-Fock method, the effect of other electrons are 
accounted for in a mean-field theory context. The orbitals are optimized by 
requiring them to minimize the energy of the respective Slater determinant. 
This operation leads to the Hartree-Fock equation described above. 
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Since the Fock operator depends on the orbitals used to construct the 
corresponding Fock matrix, the eigenfunctions of the Fock operator are in 
turn new orbitals which can be used to construct a new Fock operator. In 
this way, the Hartree-Fock orbitals are optimized iteratively until the change 
in total electronic energy falls below a predefined threshold. In this way, 
a set of self-consistent one-electron orbitals are calculated. The Hartree- 
Fock electronic wavefunction is then the Slater determinant constructed out 
of these orbitals. Following the basic postulated of quantum mechanics, the 
Hartree-Fock wavefunction can then be used to compute any desired chemical 
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or physical property within the framework of the Hartree-Fock method and 
the approximations employed. 
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Chapter 5 



An Introduction to Configuration 
Interaction Theory 

Adapted from C. D. Sherrill's notes: "An Introduction to Configuration In- 
teraction Theory". 



5.1 Introduction 

These notes attempt to present the essential ideas of configuration interac- 
tion (CI) theory in a fairly detailed mathematical framework. Of all the ab 
initio methods, CI is probably the easiest to understand, and perhaps one of 
the hardest to implement efficiently on a computer! The next two sections 
explain what the CI method is: the matrix formulation of the Schrodinger 

equation = E The remaining sections describe various simplifica- 
tions, approximations, and computational techniques. 

Much of the notation used in these notes is consistent with that of Szabo 
and Ostlund, Modern Quantum Chemistry (see references). Below are listed 
several of the commonly-used symbols and their meanings. 

N The number of electrons in the system. 

n a The number of a electrons. 

Up The number of j3 electrons. 

n The number of orbitals in the one-particle basis set. 

5ij kronecker delta function, equal to one if i — j and zero otherwise. 

H The exact nonrelativistic Hamiltonian operator. 

h The one-particle part of H , H = h i J r \ . 1 / 
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H The Hamiltonian matrix, i.e. the matrix of H, in whatever ^-electron 
basis is currently being used. 

Hij The i,j-th element of H, equal to (<&i\H\tyj), where $j and are 
iV-electron CI basis functions. 

Xj The space and spin coordinates of electron i. 

r,j The spatial coordinates of electron %. 

<pi The i-th one-particle basis function (orbital). Usually denotes a spin- 
orbital obtained from an Hartree-Fock procedure. May also be written 
simply as i. 

Xi The i-th one-particle basis function (orbital) used to expand the HF 
orbitals, 0j. Usually denotes an atom spin-orbital. 

The i-th iV-electron basis function. Usually denotes a single Slater 
determinant, but may also be a configuration state function (CSF). 

I^i) Usually denotes an eigenfunction of H. The exact nonrelativistic wave- 
function if a complete basis is used in the expansion of H 

An iV-electron basis function which differs from from some reference 
function |$ ) by the replacement of spin-orbital a by spin orbital r. 
Usually implies a single Slater determinant. 

\ab . . . c) A Slater determinant with spin-orbitals a, b, ... c occupied, 



b a <p b ... C ) 



Ux 2 ) 



Mxjv 



^ C (X 2 ) 



(i\h\j) One-electron integral in the physicists' notation (i and j) are spin- 
orbitals). More explicitly, 



(Xi)/l(xi)0 i (xi) rfXi 



One-electron integral in the physicists' notation (i and j) are spin- 
orbitals). More explicitly, 

(i\h\j) One-electron integral in the chemists' notation (i and j) are spin- 
orbitals) . 
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(ij\kl) A simple two-electron integral, in physicists' notation, where k and 
I are spin-orbitals. 

(ij\kl) = / 0*(x 1 )0*(x 2 ) 1 f fc (xi)0/(x 2 ) c/xirfx 2 

,/ |ri-r 2 | 

[z_7 1 A:Z] A simple two-electron integral, in chemists' notation, where i,j,k and 
I are spin-orbitals. 

[ij\kl] = [ 0*(xi)0 i (xi) 1 -0£(x 2 )(^(x 2 )dxidx 2 

J | ri — r 2 | 

(ijl/c/) A simple two-electron integral, in chemists' notation, where i,j, k and 
I are spatial orbitals. 

(ij\kl)= 1 0,*(ri)0 i (ri) 1 f 0^(r 2 )0/(r 2 )dricZr 2 

J | r l _ r 2| 

(ij\\kl) Asymmetrized two-electron integral, equal to (ij\kl) — (ij\lk). 



61 



5.2 Fundamental Concepts 



5.2.1 Scope of the Method 

The scope of CI is to improve the HF solution by increasing the space of 
all possible many-electron wavefunction from a single Slater determinant (in 
Hartree-Fock theory) to a set of, in principle infinite, Slater determinants. 
In the following we we learn a method for a systematic generation of Slater 
determinants starting from the one-electron orbitals obtained from the HF 
calculation 1 . 



5.2.2 What is a configuration Interaction 

We learned in chapter 2.1 that any solution of the eigenvalue equation 

Htyj(r; R) = Ej^j(r; R) (5.1) 

can be expanded as a linear combination of wavefunctions that belong to a 
complete basis set 

M 

i* i > = X) c «w ( 5 - 2 ) 

i=i 

In the case of HF theory, this expansion had a single element, the Slater 
determinant made of the occupied HF one-electron orbitals. 
More generally, an arbitrary iV-electron wavefunction can be expressed ex- 
actly as a linear combination of all possible iV-electron Slater determinants 
formed from a complete set of spin orbitals {xi(x)}. If we solve the matrix 
mechanics problem H|\P) = E\fy) in a complete basis of iV-electron functions 
as just described, we will obtain all electronic eigenstates of the system ex- 
actly. If our iV-electron basis functions are denoted the eigenvectors of 
H are given as 

M 
i 

if there are M possible iV-electron basis functions (/ will be infinite if we 
actually have a complete set of one electron x%- The matrix H is constructed 
so that Hij = (&i\H\§j) for i , j = 1,2,..., M. The matrix elements may 
be written in terms of one- and two-electron integrals according to "Slater's 
rules", as discussed in section 5.4. 

The iV-electron basis functions can be written as substitutions or 



1 Remember that the HF equation can be solved, in principle, for an infinite number of 
HF orbitals, N occupied plus an infinite number of unoccupied (virtual) orbitals 
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"excitations" from the Hartree-Fock "reference" determinant, i.e. 

\*) = c \%)+J2<\^ + E c Xi)+ E &:£)+■■■ (5-3) 

ra a<b,r<s r<s<t,a<b<c 

where means the Slater determinant formed by replacing spin-orbital a 
in |$o) with spin orbital r, etc. Every ^-electron Slater determinant can be 
described by the set of N spin orbitals from which it is formed, and this set of 
orbital occupancies is often referred to as a configuration. Thus the config- 
uration interaction method is, in its most straightforward implementation, 
nothing more than the matrix mechanics solution of the time-independent 

non-relativistic electronic Schrodinger equation = E$ . One of the great 
strengths of the CI method is its generality; the formalism applies to excited 
states, to open-shell systems, and to systems far from their equilibrium ge- 
ometries 2 . 

In practice, one does not have a complete set of one-particle basis func- 
tions (xj(x)}; typically one assumes that the incomplete one-electron basis 
set is large enough to give useful results, and the CI procedure is not modi- 
fied. The quality of the one-particle basis set can be checked by comparing 
the results of calculations using progressively larger basis sets. 

It is also possible to reduce the size of the TV-electron basis set. If we 
desire only wavefunctions of a given spin and/or spatial symmetry, as is 
usually the case, we need include only those iV-electron basis functions of 
that symmetry, since the Hamiltonian matrix is block-diagonal according to 
space and spin symmetries. If one performs the matrix mechanics calculation 
using a given set of one-particle functions {xj(x)} and all possible iV-electron 
basis functions {|$i)} (possibly symmetry-restricted), the procedure is called 
full CI. The full CI corresponds to solving Schrodinger's equation exactly 
within the space spanned by the specified one-electron basis. If the one- 
electron basis is complete (it never is in practice, but it may be in theory), 
then the procedure is called a complete CI. 

Unfortunately, even with an incomplete one-electron basis, a full CI is 
computationally intractable for any but the smallest systems, due to the 
vast number of iV-electron basis functions required. The CI space must be 
reduced somehow-hopefully in such a way that the approximate CI wave- 
function and energy are as close as possible to the exact values. The effective 
reduction of the CI space is a major concern in CI theory, and we will discuss 
some of the more popular approaches in these notes. 

By far the most common CI approximation is the truncation of the CI 
space expansion according to excitation level relative to the reference state. 
The widely-employed CI singles and doubles (CISD) wavefunction includes 
only those iV-electron basis functions which represent single or double excita- 



2 By contrast, traditional single-reference perturbation theory and coupled-cluster ap- 
proaches generally assume that the reference configuration is dominant, and they may fail 
when it is not. 
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tions relative to the reference state. Since the Hamiltonian operator includes 
only one- and two-electron terms, only singly and doubly excited configura- 
tions can interact directly with the reference, and they typically account for 
about 95% of the correlation energy in small molecules at their equilibrium 
geometries. 



5.3 The Correlation Energy 

Approximate CI methods can be judged according to what fraction of the 
correlation energy they recover. The correlation energy is defined as the 
difference between the energy in the Hartree-Fock limit (Ehf) an d the exact 
nonrelativistic energy of a system (£ ) 

-E^orr = £o — Ehf (5.4) 

This energy will always be negative because the Hartree-Fock energy is an 
upper bound to the exact energy (this is guaranteed by the variational the- 
orem). The exact nonrelativistic energy £ could, in principle, be calculated 
by performing a full CI in a complete one-electron basis set. If we have an 
incomplete one-electron basis set, then we can only compute the basis set 
correlation energy, which is the correlation energy for a given one-electron 
basis. For convenience, the basis set correlation energy is often simply re- 
ferred to as the correlation energy. 

The correlation energy is the energy recovered by fully allowing the elec- 
trons to avoid each other; Hartree-Fock improperly treats interelectron re- 
pulsions in an averaged way. However, there is some inconsistency in this 
line of thinking. When a molecule is pulled apart, the electrons shouldn't 
need to avoid each other as much, so the magnitude of the correlation energy 
should decrease. In fact, the opposite is true, as shown by the basis set corre- 
lation energies given in following Table for H2O at three different geometries. 



Geometry 


E COVT (hartree) 


R, e 


-0.148028 


1.5R e 


-0.210992 


2.0R e 


-0.310067 



The correlation energy increases at stretched geometries, because our 
definition of the correlation energy, E corr = So — Ehf, includes not only the 
concept of electrons avoiding each other, which is called the dynamical cor- 
relation energy, but also a more subtle effect called the nondynamical, static 
correlation energy. Nondynamical correlation energy reflects the inadequacy 
of a single reference in describing a given molecular state, and is due to nearly 
degenerate states or rearrangement of electrons within partially filled shells. 
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Siegbahn 3 offers the following explanation of the difference between dy- 
namical and nondynamical correlation energies: 

"In many situations it is further convenient to subdivide the correlation 
energy into two parts with different physical origins. For chemical reactions 
where bonds are broken and formed, and for most excited states, the major 
part of the correlation energy can be obtained by adding only a few extra 
configurations besides the Hartree-Fock configuration. This part of the corre- 
lation energy is due to near degeneracy between different configurations and 
has its origin quite often in artifacts of the Hartree-Fock approximation. The 
physical origin of the second part of the correlation energy is the dynamical 
correlation of the motion of the electrons and is therefore sometimes called 
the dynamical correlation energy. Since the Hamiltonian operator contains 
only one- and two-particle operators this part of the correlation energy can 
be very well described by single and double replacements from the leading, 
near degenerate, reference configurations." 



5.4 Slater's Rules 

Whether we perform a full CI or only a limited CI, we must be able to express 
H n in matrix form so that we can diagonalize it and obtain the eigenvectors 
and eigenvalues of interest. In this section we discuss Slater's rules (or the 
Slater-Condon rules), which allow us to express matrix elements 

^ = (*i\H\*i) 
in terms of one- and two-electron integrals. 

At the moment, we will express these results in terms spin-orbitals using 
physicist's notation. The one-electron integrals are written as 

(i\h\j) = j ^(r 1 )^(r 1 )^(r 1 )dr 1 (5.5) 

and the two-electron integrals are written as 

(ij\\kl) = (ij\kl) - {ij\lk) (5.6) 

(ij\kl) = [ 0*( ri U*(r 2 )— Mri)<t>ifa)d*id*2 ( 5 - 7 ) 

J r 12 

Before Slater's rules can be used, the two Slater determinants must be 
arranged in maximum coincidence. Remember that switching columns in a 



3 Pe.E. Siegbahn, The direct CI method, in Methods in Computational Molecular 
Physics, edited by G.H.F. Diercksen and S. Wilson, pages 189-207, D. Reidel, Dordrecht, 
1983. 
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determinant introduces a minus sign. For instance, to calculate ($ 1 |if|$ 2 ) 
where we have 

|$i) = \abcd) 
1^2) = \crds) 

then we must first interchange columns of or |$ 2 ) to make the two 
determinants look as much alike as possible. For example, we may rearrange 
|$2> as 

I $2) — \crds) = —\crsd) = \srcd) (5.8) 

After the determinants are in maximum coincidence, we see how many spin 
orbitals they differ by, and we then use the following rules: 

1. Identical Determinants: If the determinants are identical, then 

N N 

= y^{m\h\m) + y^(mn|[mn) (5.9) 

m m>n 

2. Determinants that Differ by One Spin Orbital: 

= \... mn ...) (5.10) 
I $2) = |---pn---) 

N 

($i|#|$ 2 ) = (m\h\p) + y^{mn\\pn) 



3. Determinants that Differ by Two Spin Orbitals: 

= \---mn---) (5.11) 
|$ 2 ) = \...pq...) 

($i\H\® 2 ) = (mn\\pq) 

4. Determinants that Differ by More than Two Spin Orbitals: 

= \---mno---) (5-12) 
| < &2) = \---pqr---) 
($i|#|$ 2 ) = 

The derivation of these rules can be found in the book of Szabo and Ostlund, 
section 2.3.4 (pp. 77-81). 
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5.5 The Solution of the CI equation: The Vari- 
ational Equation 



One particular nice feature of CI method is that the calculated lowest energy 
eigenvalue is always an upper bound to the exact ground state energy. 
This is a direct consequence of the variational theorem (see sections 2.3.3 
and 2.3.4) applied to the total energy expectation value, 

E ^^ff (5.13) 

where \^>q) is an approximate solution given as a linear combination of N- 
electron wavefunction, 

l*o> = c |$o> + $>K>+ £ csi*2)+ + ■■■ 

ra a<b,r<s r<s<t,a<b<c 

or, using a more compact notation, 

|*o) = £c 04 |$ l ) (5.14) 

8=0 



Once the expansion basis is chosen, {|$j)}, the resultion secular equation 
(see sections 2.3.3 and 2.3.4), 

HC = EC (5.15) 

can be solved numerically. The outcome of this calculation is the ground 
state energy and the coefficients, {coi} for the expansion of the ground 
state iV-electron wavefunction. Note that quality of the results depend on 
the quality and the number number of the chosen basis {| $,;)}. 

At this point it is reasonable to ask why we wish to minimize the energy 
by varying the coefficients in equation 5.14. How do we know that this will 
give us the best estimate of the wavefunction? There are two answers to 
this. First, as we have just shown, minimizing the energy by variation of the 
linear expansion coefficients gives the Schrodinger equation in matrix form; 
thus the procedure is justified a posteriori by the validity of its result. 
The other reason is that, for the ground state, the linear expansion in equa- 
tion gives an expectation value for the energy E, which is always an upper 
bound to the exact ground state energy (this is called the Variational The- 
orem). The best estimate of E, then, is the minimum value which can be 
obtained by varying the coefficients in equation 5.14 (while also maintaining 
normalization). These arguments also hold for excited states, so long as each 
excited state is made orthogonal to all lower states. 
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Proof of the Variational Theorem for the Ground State. 
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5.6 Classification of Basis Functions by Excita- 
tion Level 



Now we will discuss the importance of various excitation classes to the CI 
wavefunction. As noted in the equation 

l*> = col$o> + $>^>+ £ C X) + E <££!*&> + ■■■ 

ra a<b,r<s r<s<t,a<b<c 

the CI expansion is typically truncated according to excitation level; in the 
vast majority of CI studies, the expansion is truncated (for computational 
tractability) at doubly-excited configurations. Since the Hamiltonian con- 
tains only two-body terms, only singles and doubles can interact directly 
with the reference. This is a direct result of Slater's Rules. The structure of 
the CI matrix with respect to excitation level is given below (adapted from 
Szabo and Ostlund, p. 235), where \S), \D), \T), represent blocks of singly, 
doubly, triply, and quadruply excited determinants, respectively. The Hamil- 
tonian matrix H is Hermitian; if only real orbitals are used, as is usually the 
case, then the Hamiltonian is also symmetric. Thus only the lower triangle 
of H is shown below. 



o 



<<&, 
(S\ 

H = T 

(Q\ 



(®o\H\$ ) 

(S\H\S) 

(D\H\$ ) (D\H\S) (D\H\D) 

o \t\h\s) \t\h\d) (T\H\T) 

(Q\H\D) (Q\H\T) {Q\H\Q} 



Note that the matrix elements (S\H\Qq) are given as 0. This is due to Bril- 
louin's theorem, which is valid when our reference function |$o) is obtained 
by the Hartree-Fock method (Hartree-Fock guarantees that off-diagonal ele- 
ments of the Fock matrix are zero, and it turns out that the matrix element 
between two Slater determinants which differ by one spin orbital is equal 
to an off-diagonal element of the Fock matrix). Furthermore, the blocks 
(X\H\Y) which are not necessarily zero may still be sparse; for example, the 
matrix element ) which belongs to the block (D\H\Q), will be 

nonzero only if a and b are contained in the set {c, d, e, /} and if r and s are 
contained in the set {t, u, v, w}. 

Since only the doubles interact directly with the Hartree-Fock reference, 
we expect double excitations to make the largest contributions to the CI 
wavefunction, after the reference state. Indeed, this is what is observed. 
Even though singles, triples, etc. do not interact directly with the refer- 
ence, they can still become part of the CI wavefunction (i.e. have non-zero 
coefficients) because they mix with the doubles, directly or indirectly. Al- 
though singles are much less important to the energy than doubles, they are 
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generally included in CI treatments because of their relatively small number 
and because of their greater importance in describing one-electron properties 
(dipole moment, etc.) 



5.7 Energy Contributions of the Various Exci- 
tation Levels 

The following table demonstrates the importance of various excitation classes 
in obtaining CI energies. We see that singles and doubles account for 95% of 
the correlation energy at the equilibrium geometries of the molecules listed. 
We see that quadruple excitations are more important than triples, at least 
as far as the energy is concerned. At stretched geometries, the CISD and 
CISDT methods become markedly poorer, yet the CISDTQ method still re- 
covers a very high (and nearly constant) fraction of the correlation energy, 
suggesting that CISDTQ should give reliable results for energy differences 
across potential energy surfaces for molecules of this size. 



Percent Correlation Energy 



Molecule 


CISD 


CISDT 


CISDTQ 


BH 


94.91 


n/a 


99.97 


H 2 0(R e ) 


94.70 


95.47 


99.82 


H 2 0(1.5R e ) 


89.39 


91.15 


99.48 


H 2 O(2.0R e ) 


80.51 


83.96 


98.60 


NH 3 


94.44 


95.43 


99.84 


HF 


95.41 


96.49 


99.86 




96.36 


96.87 


99.96 



5.8 Size of the CI Space as a Function of Exci- 
tation Level 

From the next table, we can see that the number of iV-electron basis func- 
tions increases dramatically with increasing excitation level. It should be 
pointed out that while the calculations on BH, HF, and hy used DZP basis 
sets, those on h 2 and NH 3 used only DZ basis sets. A DZP basis should 
be considred the minimum adequate basis for a truly meaningful benchmark 
study, and even then it is desirable to use a high-quality basis such as an 
Atomic Natural Orbital (ANO) set. 



70 



CFS's required 


Molecule 


C1SD 


C1SDT 


C1SDTQ 


FC1 


BH 


568 


n/a 


28'698 


132'686 


H 2 


361 


3'203 


17'679 


256'473 


NH 3 


461 


4'029 


19'925 


137'321 


HF 


552 


6712 


48'963 


944'348 


H 7 + 


V271 


24'468 


248'149 


2'993'933 



While it is generally possible to perform CISD calculations on small 
molecules with a good one-electron basis, the CISDTQ method is limited 
to molecules containing very few heavy atoms, due to the extreme number of 
iV-electron basis functions required. Full CI calculations are of course even 
more difficult to perform, so that despite their importance as benchmarks, 
few full CI energies using flexible one-electron basis sets have been obtained. 

The size of the full CI space in CSF's can be calculated (including spin 
symmetry but ignoring spatial symmetry) by Weyl's dimension formula as 
applied to the Distinct row table (DRT). If N is the number of electrons, n 
is the number of orbitals, and S is the total spin, then the dimension of the 
CI space in CSF's is given by 

2S + 1 ( n + l \ ( n + 1 \ 
VnNS ~ \N/2-S ){N/2 + S + l J 

The dimension of the full CI space in determinants (again, ignoring spatial 
symmetry) is computed simply by 



or, equivalently, 



n \ I n 

N a [ N 



D 



nNS 



n \ n 
N/2 + S J I N/2 - S 
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Chapter 6 

Many-Body Perturbation Theory 



The idea in perturbation methods is that the problem under investigation 
only differs slightly from a problem which has already been solved (exactly 
or approximately). The solution to the given problem should therefore in 
some sense be close to the solution of the already known system. This is 
described mathematically by defining a Hamiltonian operator which consists 
of two parts, a reference (Ti-o) and a perturbation (T~C). 

The basic assumption of perturbation methods is that the 7i' operator in 
some sense is "small" compared to H,q. In quantum mechanics, perturbation 
methods can be used for adding corrections to solutions which employ an 
independent particle approximation, and the theoretical framework is then 
called Many-Body Perturbation Theory (MBPT). 



6.1 Perturbation Theory in Quantum Mechan- 
ics 

Let us assume that the Schrodinger equation for the reference Hamiltonian 
operator is 

7i = H-o + A Tt 

Wo$* = i = 0,1, 2,3,..., oo (6.1) 

The solutions for the unperturbed Hamiltonian operator form a complete 
set (since 7i is hermitian) which can be chosen to be orthonormal, and 
A is a (variable) parameter determining the strength of the perturbation. 
The parameter A can be varied systematically to switch the system from the 
unperturbed case (A = 0) to the fully perturbed (A = 1) case. At present we 
will only consider cases where the perturbation is time-independent, and the 
reference wave function is nondegenerate. To keep the notation simple, we 
will furthermore only consider the lowest energy state (i = 0). 
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The perturbed Schrodinger equation is 

Wtf(A) =E(X)%(\) (6.2) 

If A = 0, then H = H , *(0) = $ and E(0) = E°. As the perturbation is 
increased from zero to a finite value, the new energy and wavefunction must 
also change continuously, and they can be written as a Taylor expansion in 
powers of the perturbation parameter A (from now on, we will consider the 
perturbation of the ground state orbital only, $ — 

E(X) = A E<® + A 1 ^ + X 2 E^ + X 3 E^ + ... 

tf(A) = A * (0) + A^ (1) + A 2 * (2) + A 3 * (3) + . . . (6.3) 

For A = 0, it is seen that A = 1, * (0) = $ and E^ = E°, these are 
the unperturbed, or zero-order wavefunction and energy. The ^f^\ ^^ 2 \ . . . 
and E^,E^ 2 \... are the first-, second-, etc. order corrections. The A 
parameter will eventually be set equal to 1, and the nth order energy or 
wavefunction become a sum of all terms up to order n. 



6.1.1 Normalization condition 

It is convenient to choose the perturbed wavefunction so that the overlap 
with the unperturbed wavefunction is equal 1. This has the consequence 
that all correction terms are orthogonal to the reference wavefunction. We 
therefore have 

($|$) = 1 

+ A^ (1) + AV 2) ... |$) = 1 (6.4) 
+ A 1 ^ + A 2 (* (2) |$) + • • • = 1 

which implies 

= o (6.5) 



6.1.2 The nth-order perturbation equation 

Once all the correction terms have been calculated, it is trivial to normalize 
the total wavefunction. 

With the expansions (Eq.6.3) the Schrodinger equation (6.2) becomes 
(Ho + XH')(X°¥ 0) + A^ (1) + X 2 ¥ 2) + ...) = 

( A °£(°) + A 1 ^ 1 ' + X 2 E^ + . . . ) (A V°> + A 1 ^ + A V 2 > + . . . ) (6.6) 
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Since this holds for any value of A, we can collect terms with the same power 
of A to give 

A (0) . Ho ^(0) = £(0)^(0) 

A« : H + OT (0) = + 

A< 2 >: H ^ (2) + = £ (0) V® + + £ (2) *(°) (6.7) 

A (n) . Hq ^(n) + = ^ E (j) y(n-j) 

3=1 

These are zero-, first-, second-, nth-order perturbation equations. The zero- 
order equation is just the Schrodinger equation for the unperturbed problem. 
The first-order equation contains two unknowns, the first-order correction to 
the energy, E^\ and the first-order correction to the wavefunction, 



6.1.3 Ray leigh- Schrodinger perturbation formula 

Up to this point we are still dealing with undetermined quantities, energy and 
wavefunction corrections at each order. The first-order equation is one equa- 
tion with two unknowns. Since the solutions to the unperturbed Schrodinger 
equation generates a complete set of orthogonal functions, the unknown first- 
order correction to the wave function, can be expanded in these func- 
tions, {$i}°^ - This is known as Rayleigh-Schrodinger perturbation 
theory. 



First order perturbation. The A 1 equation in (6.7) becomes 

(H - E<®) + (W - £ (1) )$ = (6.8) 

with 

oo 

#W = ^ Ci ^ ( 6 .9) 

and yields (exercise) 

= ($ |ft'|$ ) (6.10) 

This shows that the first-order correction to the energy is an average of the 
perturbation operator over the unperturbed wave function. 
The first-order correction to the wavefunction can be obtained by multiplying 
Eq. 6.7 for A 1 from the left by the function corresponding to a given Cj (<&j) 
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and integrating to give 

The expansion coefficients determine the first-order correction to the per- 
turbed wavefunction (Eq.6.8), and they can be calculated using the known 
unperturbed wavefunctions and energies. The coefficient in front of $ f° r 
(Eq.6.9) cannot be determined from the above formula, but the assump- 
tion of intermediate normalization (Eq.6.4) makes c = 0. 



Second order perturbation. Starting from the second-order perturba- 
tion equation (Eq. 6.7) and using intermediate normalization (co = do = 0), 
the second order energy correction is expressed as 

vD(2) = ^ di $. ( 6 .i2) 

i 

with 

^^^^ MfflMW (6 . 13) 

and second-order wavefunction expansion coefficients 

J ^ (E -^)(i?o-^) (^o-^) 2 1 ' ' 

The formulas for higher-order corrections become increasingly complex. The 
main point, however, is that all corrections can be expressed in terms of 
matrix elements of the perturbation operator over the unperturbed wave- 
functions, and the unperturbed energies. 



6.2 M0ller-Plesset Perturbation Theory 

So far the theory has been completely general. In order to apply perturbation 
theory to the calculation of correlation energy, the unperturbed Hamiltonian 
operator must be selected. 

The most common choice is to take the sum of the single particle Fock op- 
erators as the upperturbed Hamiltonian, and the difference between the full 
Hamiltonian and the HF Hamiltonian as perturbation. This is called M0ller- 
Plesset (MP) perturbation theory. 

The sum of Fock operators counts the (average) electron-electron repul- 
sion twice (Eqs. 4.7) and the perturbation becomes the exact r^ 1 = Vij 



76 



operator minus twice the vf F operator, 



vg* =J lJ -K ij (6.15) 
= ^2<ik\jk) - (ik\kj) = E<**» ■ (6.16) 



The operator associated with this difference is often referred to as the fluc- 
tuation potential. 



The starting point is the Hartree-Fock operator, Eq. 4.20, (i and j are 
electron indices) 

N 

^ = ^ + ^(J ? -/C J ) (6.17) 

j 

with 

Fi<t>i{xi) = £iM x i) (6-18) 
The Coulomb and the exchange operators are defined as 

Jj\(pi(xi)) = (0j(x)|uy|^-(a:))|0i(xi)) (6.19) 
K.j\(f)i(xi)) = {(j)j{x)\vij\(j)i{x))\(j)j{xi)) (6.20) 

Starting from Ti (the operator for the electron with index i) we can construct 
the "N-electrons" operator Yli=i Fii w hich is taken as reference Hamiltonian 

N N N N N N 

«o = ££ = 5>+££ (4 - *«) = E + E' i! S F ■ ( 6 - 21 ) 

i=l i=l i=l j'=l i=l M=l 

The perturbation Hamiltonian is therefore given by 

N N 

H' = H-H = Y;^-Yltf F (6.22) 
where for a given Slater determinant many-electron wavefunction $o ex_ 
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pressed in the HF orbital basis {0j} (compare with Eq. 4.7 1 ) 

N N 

<$o|$>^°> = 2< $ °lE'C>o) = (Vee) (6.23) 

The zero-order wave function is the HF determinant, and the zero-order 
energy is just a sum of MO energies. 

N N 

e^ = j2(^m = J2^ f ( 6 - 24 ) 

i i 

The first-order energy correction is the average of the perturbation op- 
erator over the zero-order wave function (Eq.6.10): 



N N 

E {1) = (<W1$o> = (<&olX>il$o) " (*o\J^ v* F \*o) (6-25) 

i<j i,j=l 

= (V ee )-2(V ee ) = -(V ee ) 

This yields a correction for the double counting of the electron-electron re- 
pulsion at zero order. Comparing Eq.6.22 with the expression for the total 
energy in Eq.4.7, it is seen that the first-order energy (sum of E^ and E^) 
is exactly the HF energy. 

Using the notation E(MPn) to indicate the correction at order n, and MPn 



1 For the set of HF orbitals, {<fii}, we have 

(Wo) = = E< l N*> + £»1 - MH] 

i i ij 

i ij 

i " ij 

i ij 

= J2(i\h\i) + -J2j u -K u 
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to indicate the total energy up to order n, we have 

N 

MPO : £(MP0)=^^ F 

a 

MP1 : £(MP0) + £7(MP1) = E(HF) 

Electron correlation energy thus starts at order 2 with this choice 
of H . 

In developing perturbation theory it was assumed that the solutions to 
the unperturbed problem formed a complete set. This in general means that 
there must be an infinite number of functions, which is impossible in ac- 
tual calculations. The lowest energy solution to the unperturbed problem is 
the HF wave function, additional higher energy solutions are excited Slater 
determinants, analogously to the CI method. When a finite basis set is em- 
ployed it is only possible to generate a finite number of excited determinants. 
The expansion of the many-electron wave function is therefore truncated. 



The second-order correction to the energy, which is the first contribution 
to the correlation energy, involves a sum over doubly excited determinants 
(singly excited Slater determinants give no contribution to the energy. This 
is known as Brillouins theorem). These can be generated by promoting two 
electrons from occupied orbitals a and b to virtual orbitals r and s. 
The summation must be restricted so that each excited state is only counted 
once 

£(2 )_^ ($ol^lJJ(gHl$o) (626) 



a<b r<s 



The matrix elements between the HF and a doubly excited state are given 
by two electron integrals over MOs. The difference in total energy between 
two Slater determinants becomes a difference in MO energies (essentially 
Koopmans' theorem), and the explicit formula for the second-order M0ller- 
Plesset correction is 

£(MP2) = J2 j£ ^"^^'^ - (^H^-)] 2 ( g 27 ) 

Once the two-electron integrals over MOs are available, the second-order en- 
ergy correction can be calculated as a sum over such integrals. There are of 
the order of M 4 integrals, thus the calculation of the energy (only) increases 
as M 4 with the system size. However, the transformation of the integrals 
from the AO to the MO basis grows as M 5 . MP2 is an M 5 method, but 
fairly inexpensive as not all two-electron integrals over MOs are required. 
Only those corresponding to the combination of two occupied and two vir- 
tual MOs are needed. In practical calculations this means that the MP2 
energy for systems with 100-150 basis functions can be calculated at a cost 
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similar to or less than what is required for calculating the HF energy. MP2 
typically accounts for 80-90% of the correlation energy, and it is the most 
economical method for including electron correlation. 



Some remarks 

The main limitation of perturbation methods is the assumption that the zero- 
order wave function is a reasonable approximation to the real wave function, 
i.e. the perturbation operator is sufficiently "small". The poorer the HF 
wave function describes the system, the larger are the correction terms, and 
more terms must be included to achieve a given level of accuracy. If the 
reference state is a poor description of the system, the convergence may be 
so slow or erratic that perturbation methods cannot be used. 

Actually it is difficult to prove that the perturbation expansion is conver- 
gent, although many systems show a behavior which suggests that it is the 
case. This may to some extent be deceptive, as it has been demonstrated that 
the convergence properties depend on the size of the basis set and the major- 
ity of studies have employed small or medium sized basis sets. A convergent 
series in a DZP type basis for example may become divergent or oscillating 
in a larger basis, especially if diffuse functions are present. In the ideal case 
the HF, MP2, MP3 and MP4 results show a monotonic convergence towards 
a limiting value, with the corrections being of the same sign and numerically 
smaller as the order of perturbation increases. Unfortunately, this is not the 
typical behavior. Even in systems where the reference is well described by 
a single determinant, oscillations in a given property as a function of per- 
turbation order are often observed. This is not completely understood, but 
may at least partly be due to the fact that the choice of the unperturbed 
Hamilton operator does not make the perturbation particularly small. 

In practice only low orders of perturbation theory can be carried out, 
and it is often observed that the HF and MP2 results differ considerably, the 
MP3 result moves back towards the HF and the MP4 moves away again. For 
"well-behaved" systems the correct answer is normally somewhere between 
the MP3 and MP 4 results. MP2 typically overshoots the correlation effect, 
but often gives a better answer than MP3, at least if medium sized basis sets 
are used. Just as the first term involving doubles (MP2) tends to overestimate 
the correlation effect, it is often observed that MP4 overestimates the effect 
of the singles and triples contributions, since they enter the series for the first 
time at fourth order. 
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Chapter 7 
Coupled Cluster 



Coupled cluster (CC) is nowadays one of the most prevalent methods in 
quantum chemistry that includes electronic correlation. 

Perturbation methods (like MPn) add all type of corrections (S,D,T, 
etc) to the reference wavefunction up to a given order (2,3,4, etc) 1 . 
The idea in Coupled Cluster (CC) methods is to include all corrections of 
a given type (S,D,T, etc) to infinite order. 

The wavefunction and the ground state energy are denoted by and E, 
respectively. Other variants of the coupled-cluster theory, such as equation- 
of-motion coupled cluster and multi-reference coupled cluster may also pro- 
duce approximate solutions for the excited states of the system. 

Properties: 

• CC in non-variational 

• CC is size-extensive 

• CC is size-consistent if the reference wavefunction is size consistent 

The wavefunction of the coupled-cluster theory is written as an exponen- 
tial ansatz: 

|^) = e t |$ ) (7.1) 
where |$o) is a Slater determinant usually constructed from Hartree-Fock 

molecular orbitals. The so called cluster operator, T, is an excitation op- 
erator which, when acting on \&o), produces a linear combination of excited 
Slater determinants. 

The choice of the exponential ansatz is opportune because (unlike other 



The order of a correction corresponds to the exponent in the perturbation expansion. 
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ansatzes, for example, configuration interaction) it guarantees the size con- 
sistency 2 and size extensivity 3 of the solution. 

7.1 The cluster operator 

The cluster operator is written in the form 

f = f 1 + f 2 + f 3 + ... (7.2) 

where T\ is the operator of all single excitations, T 2 is the operator of all 
double excitations and so forth. In the formalism of second quantization 
these excitation operators are conveniently expressed as 

a r 

T2 = ^ tr ab^ a a b alal (7.4) 

a,b r,s 

In the above formulae a* and a denote the creation and annihilation opera- 
tors respectively and a, b stand for occupied and r, s for unoccupied orbitals. 
The creation and annihilation operators in the coupled cluster terms above 
are written in canonical form, where each term is in normal order (i.e. with 
all creation operators on the right hand side). The one-particle excitation 

operator and the two-particle excitation operator, T\ and T 2 convert the refer- 
ence function |$o) m t° a linear combination of the singly- and doubly-excited 
Slater determinants, respectively. Solving for the unknown coefficients t r a and 
t r a l is necessary for finding the approximate solution | . 

Taking into consideration the structure of T, the exponential operator e T 
may be expanded into a Taylor series: 

rp2 Ji3 Ji2 rp2 

e T = l+f + — + — + ... = l+f 1 +f 2 + ^-+f 1 f 2 + ^- + ... (7.5) 
This series is finite in practice because the number of molecular orbitals is 



2 Size consistency (or strict separability) is a property that guarantees the consistency 
of the energy behavior when interaction between the involved molecular system is nullified 
(for example, by increasing the distance between two moieties: E(AB) has the right limit, 
E(A) + E(B) for large separation of the two molecules, A and B.) 

3 In physics and chemistry an intensive property (also called a bulk property) of a 
system is a physical property of the system that does not depend on the system size or the 
amount of material in the system. By contrast, an extensive property of a system does 
depend on the system size or the amount of material in the system. 
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finite, as is the number of excitations. In order to simplify the task for finding 
the coefficients t^'", the expansion of into individual excitation operators is 

terminated at the second or slightly higher level of excitation (rarely exceed- 
ing four). This approach is warranted by the fact that even if the system 

admits more than four excitations, the contribution of T 5 , T 6 etc to the oper- 
ator T is small. Furthermore, if the highest excitation level in the T operator 
is n, 

f = 1 + Ti + • ■ ■ + f n (7.6) 

then Slater determinants excited more than n times may (and usually do) 
still contribute to the wave function because of the non-linear nature of 

the exponential ansatz (i.e. for n = 2 the expansion contains terms like TiT 2 

and Tf , which correspond to triple and quadruple excitations, respectively). 

Therefore, coupled cluster terminated at T n usually recovers more correlation 
energy than configuration interaction with maximum n excitations. 



7.2 The coupled cluster energy 



With the coupled cluster wave function the Schrodinger equation becomes 

He f |$ ) = Ee f |$ ) (7.7) 
Multiplying from the left by (<3>o| and integrating gives 

(<S> \He f |$ ) = £ cc ($o|/$ > (7.8) 

= £ cc ($ |(i+Ti+T 2 + ...)$o> (7-9) 

and therefore 

E cc = (%\He f |$ ) (7.10) 
(all terms ($ |Tj|$ ) are zero. Why?). 

In the case in which excitations are limited to single and double excited Slater 
determinants, the cc energy simplifies to 

occ vir occ vir 

a r ab(a<b) rs(r<s) 

(7.11) 

When using HF orbitals for constructing the Slater determinants, the first 
matrix elements are zero (Brillouins theorem) and the second matrix elements 
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are just two-electron integrals over molecular orbitals, 

occ vir 

E CC = E + E (C6+*a*6-OD((0a0 6 |^|0r0 S )-(0a0 6 |^|0 S 0r)) (7.12) 

ab(a<b) rs(r<s) 

What is still missing are the expansion amplitudes, t™ b \. 

7.3 The coupled cluster equations for ampli- 
tudes, t r a l"' 

Coupled-cluster equations are equations whose solution is the set of coeffi- 
cients t r a l"\ There are several ways of writing such equations but the standard 
formalism results in a closed set of equations which may be solved iteratively. 
The interested students are invited to read the literature (see for example 
the book of Szabo and Ostlund, pp. 286-319). 

7.4 Types of coupled-cluster methods 

The classification of traditional coupled-cluster methods rests on the highest 

number of excitations allowed in the definition of T. The abbreviations for 
coupled-cluster methods usually begin with the letters "CC" (for coupled 
cluster) followed by 

S - for single excitations (shortened to singles in coupled-cluster termi- 
nology) 

D - for double excitations (doubles) 

T - for triple excitations (triples) 

Q - for quadruple excitations (quadruples) 

Thus, the T operator in CCSDT has the form 

f = f i + X 2 + f ' 3 

Terms in round brackets indicate that these terms are calculated based on 
perturbation theory. For example, a CCSD(T) approach simply means: 

1. A coupled-cluster method 

2. It includes singles and doubles fully 

3. Triples are calculated with perturbation theory. 
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Chapter 8 

Density functional theory 



The content of this chapter is partially adapted from: 

Klaus Capelle, A Bird's-Eye View of Density-Functional Theory, arXivxond- 
mat/0211443v5 

The electron density is a more attractive quantity 1 , than many particle 
wave function which depends on all coordinates of all particles, i.e., for N 
electrons, it depends on 3N variables (or 4N if you count in spin). 

<P(xi,x 2 , . . . ,Xjv) = M fc det |^a(xi)^ b (x 2 ) . . .i[)x(xn)\ 



p(xi) = M J ... / ^(xi, x 2 , . . . , xat) ^*(xi, x 2 , . . . , xat) dx 2 . . . dx N 
where M is a normalization factor. 

Why DFT? 

A practical reason 

The complete oxygen wavefunction, ^(xi,x 2 , . . . , Xjv), (8 electrons) depends on 24 
coordinates. To store this wavefunction we would need 

10 entries for each coordinate > 10 24 entries 

8 byte per entry ► 10 25 bytes 

5 x 10 9 bytes per DVD ► 2 x 10 14 DVDs 

10 g per DVD ► 2 x 10 9 t DVDs 

1 depends only on x, y, z, and eventually, there may be two densities for spin polarized 
systems, one for spin up electrons p-f(r) and one for spin down electrons Pi(r) 
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Instead, in DFT we will need about 8 x 10 3 bytes to store the density of the system, 
which requires less than a floppy disk! 

A physical reason 

DFT has many advantages compared to other ab initio techniques 

• DFT is computationally very efficient. 

• DFT is conceptually simple. 

"A great strength of DFT language is its appropriateness for defining and 
elucidating important universal concepts of molecular structure and molec- 
ular reactivity. In traditional quantum chemistry this has, of course, also 
been a major goal, but it is tortuous to try to conceptualize how many-body 
wavefunctions are related to structure and behavior. In DFT not only is the 
electron density itself very easy to visualize, but there is the big advantage 
that the electron number N has a central place in the theory. After all, much 
of the chemistry is about transfer of electrons from one place to another." 
(W. Kohn,A.D.Becke, R.G. Parr, J. Phys. Chem, 12974-12980,1996.) 
Concepts like chemical potential (electronegativity) , hardness, polarizability, 
response functions and reactivity functions (Fukui functions) emerge natu- 
rally from the DFT energy functional, E[p], and its derivatives, and their 
computation becomes possible. 

• DFT can be easily combined with molecular dynamics of the nuclei. 
Forces on classical ions can be computed using the Hellmann-Feynman the- 
orem starting from the electronic density. 



8.1 What is density- functional theory 

Density-functional theory is one of the most popular and successful quan- 
tum mechanical approaches to matter. It is nowadays routinely applied for 
calculating, e.g., the binding energy of molecules in chemistry and the band 
structure of solids in physics. Also applications relevant for fields tradition- 
ally considered more distant from quantum mechanics, such as biology and 
mineralogy are beginning to appear. Superconductivity, atoms in the focus 
of strong laser pulses, relativistic effects in heavy elements and in atomic nu- 
clei, classical liquids, and magnetic properties of alloys have all been studied 
with DFT. 

In quantum mechanics we learn that all information we can possibly have 
about a given system is contained in the system's wave function, which is 
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the lowest energy solution of the many-electron Schrodinger equation 



or, in a more compact notation, 

[f e (r) + t> eJV (r, R) + VWat(R) + t> ee (r)]*(r, R) = £ e ^(r, R) (8.2) 

Convention. In the following we use v(r, R) for the external potential 
generated by the nuclei at position R/ (R = (R 1; R 2 , . . . , Rat m )- Since we 
usually consider the nuclei kept fixed in space, we will sometimes drop the 
dependence on R, and write simply v(r) = v(r, R). 

The usual quantum-mechanical approach to the Schrodinger equation 
(SE) can be summarized by the following sequence 

v(r,R) => *(r i, r2, . . . , r n ) ==k> observables 

i.e., one specifies the system by choosing v (r, R), plugs it into the Schrodinger 
equation, solves that equation for the wave function and then calculates 
observables by taking expectation values of operators with this wave function. 
One among the observables that are calculated in this way is the electron 
density, p(r). 

In the previous chapters, we have investigated methods for the solution 
of the many-electron Schrodinger equation using a systematic expansion in 
Slater determinants. The problem with these methods is the great demand 
they place on one's computational resources: it is simply impossible to apply 
them efficiently to large and complex systems. Nobody has ever calculated 
the chemical properties of a 100-atom molecule with full CI 2 . 

DFT provides a viable and accurate 3 alternative to post-HF methods 



2 A simple estimate of the computational complexity of this task is to imagine a real 
space representation of $ on a mesh, in which each coordinate is discretized by using 20 
mesh points (which is not very much). For TV electrons, ^ becomes a function of 3N 
coordinates (ignoring spin, and taking "J to be real), and 2Q 3N values are required to 
describe ^ on the mesh. The density n(r) is a function of three coordinates, and requires 
20 3 values on the same mesh. CI and the Kolm-Sham formulation of DFT additionally 
employ sets of single-particle orbitals. TV such orbitals, used to build the density, require 
20 3 N values on the same mesh. (A CI calculation employs also unoccupied orbitals, and 
requires more values.) For N — 10 electrons, the many-body wave function thus requires 
20 30 / 20 3 w 10 35 times more storage space than the density, and 20 30 /(10 x 20 3 ) « 10 34 
times more than a set of single-particle orbitals. Clever use of symmetries can reduce these 
ratios, but the full many-body wave function remains unaccessible for real systems with 
more than a few electrons. 

3 Accuracy is a relative term. As a theory, DFT is formally exact. Its performance 



= £ei*(r,R) 
(8.1) 
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for the computation of medium (100 electrons) to big size (thousands of 
electrons) systems. 

DFT explicitly recognizes that nonrelativistic Coulomb systems differ 
only by their potential v(r, R), and supplies a prescription for dealing with 

the universal operators T e and V ee once and for all. Furthermore, DFT pro- 
vides a way to systematically map the manybody problem, with V ee , onto 

a single-body problem, without V ee . All this is done by promoting the par- 
ticle density p(r) from just one among many observables to the status of 
key variable, on which the calculation of all other observables can be based. 
This approach forms the basis of the large majority of electronic-structure 
calculations in physics and chemistry. Much of what we know about the elec- 
trical, magnetic, and structural properties of materials has been calculated 
using DFT, and the extent to which DFT has contributed to the science 
of molecules is reflected by the 1998 Nobel Prize in Chemistry, which was 
awarded to Walter Kohn, the founding father of DFT, and John Pople, who 
was instrumental in implementing DFT in computational chemistry. 
The density-functional approach can be summarized by the sequence 

p(r) =^( ri ,r 2 ,...,r n ) =^(r) (8.3) 

i.e., knowledge of p(r) implies knowledge of the wave function and the poten- 
tial, and hence of all other observables. This is valid for a given fixed position 
of the nuclei and therefore we remove the explicite dependence on R from 
all quantities). Although this sequence describes the conceptual structure of 
DFT, it does not really represent what is done in actual applications of it, 
which typically proceed along rather different lines, and do not make explicit 
use of many-body wave functions. 

The literature on DFT is large, and rich in excellent reviews and overviews. 
Some representative examples of full reviews and systematic collections of re- 
search papers are: 

1. R. M. Dreizler and E. K. U. Gross, Density Functional Theory, Springer, 
Berlin, 1990. 

2. R. G. Parr and W. Yang, Density- Functional Theory of Atoms and 
Molecules, Oxford University Press, Oxford, 1989. 

3. W. Koch and M. C. Holthausen, A Chemist's Guide to Density Func- 
tional Theory, John Wiley & Sons, New York, 2001. 



in actual applications depends on the quality of the approximate density functionals em- 
ployed. For small numbers of particles, or systems with special symmetries, essentially 
exact solutions of Schrodinger equation can be obtained, and no approximate functional 
can compete with exact solutions. For more realistic systems, modern (2007) sophisti- 
cated density functionals attain rather high accuracy. Bond-lengths of molecules can be 
predicted with an average error of less than 0.01, lattice constants of solids with an average 
error of less than 0.05, and molecular energies to within less than 5 kcal/mol. (For com- 
parison: already a small molecule, such as water, has a total energy of 48'OOOfcca^/mc^). 
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4. R. O. Jones and 0. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989). 

5. J. M. Seminario (Ed.), Recent Developments and Applications of Mod- 
ern DFT, Elsevier, Amsterdam, 1996. 



8.2 Functionals and their derivatives 

Before we discuss density-functional theory more carefully, let us introduce 
a useful mathematical tool. Since according to the above sequence the wave 
function is determined by the density, we can write \1/ = \l/[p](ri, r 2 , . . . , r^) 
which indicates that \l/is a function of its N spatial variables, but a functional 
of p(r). A functional F[p] can be defined (in an mathematically sloppy way) 
as a rule for going from a function to a number, just as a function y = f(x) 
is a rule (/) for going from a number (x) to another number (y). A simple 
example of a functional is the particle number, 



which is a rule for obtaining the number N, given the function p(r). Note 
that the name given to the argument of p is completely irrelevant, since the 
functional depends on the function itself, not on its variable. Hence we do 
not need to distinguish F[p(r)] from, e.g., F[p(r')}. Another important case 
is that in which the functional depends on a parameter, such as in 



which is a rule that for any value of the parameter r associates a value 
v#[p](r) with the function p(r'). This term is the so-called Hartree potential, 
which we have already encountered in the HF theory. 

8.2.1 Functional variation 

Given a function of one variable, y = f(x), one can think of two types of 
variations of y, one associated with x, the other with /. For a fixed functional 
dependence f(x), the ordinary differential dy measures how y changes as a 
result of a variation x — > x+dx of the variable x. This is the variation studied 
in ordinary calculus. Similarly, for a fixed point x, the functional variation 
5y measures how the value y at this point changes as a result of a variation in 
the functional form f(x). This is the variation studied in variational calculus. 

8.2.2 Functional derivative 

The derivative formed in terms of the ordinary differential, df/dx, measures 
the first-order change of y = f(x) upon changes of x, i.e., the slope of the 




(8.4) 




(8.5) 
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function f(x) at x: 



fix + dx) = fix) + ^-dx + Oidx 2 ) (8.6) 
dx 

The functional derivative measures, similarly, the first-order change in a func- 
tional upon a functional variation of its argument: 

F[f(x) + 8f(x)} = F[f(x)} + J dxs(x)5f(x)dx + 0(5 f) , (8.7) 

where the integral arises because the variation in the functional F is deter- 
mined by variations in the function at all points in space. The first-order 
coefficient (or 'functional slope') s(x) is defined to be the functional deriva- 
tive SF[f]/6f(x), 

•W - (^) 

The functional derivative allows us to study how a functional changes upon 
changes in the form of the function it depends on. Detailed rules for calcu- 
lating functional derivatives are described in Appendix A of the book of Parr 
and Yang. 

A general expression for obtaining functional derivatives with respect to 
p(r) of a functional F[p] = J drf(p, p', p", p'", . . . ; r), where primes indicate 
ordinary derivatives of p(r) with respect to r, is 

SF[P] = d£_ ±d£ + £_df_ _ _^9f_ + (g g) 

5p(x) dp dx dp' dx 2 dp" dx 3 dp'" 

This expression is frequently used in DFT to obtain exchange-correlation 
(xc) potentials from xc energies. 



8.3 The Hohenberg-Kohn theorem 

At the heart of DFT is the Hohenberg-Kohn (HK) theorem. This theorem 
states that for ground states the equation 

p(r) = M J ... J ^(r, r 2 , . . . , v N ) #*(r, r 2 , . . . , r N ) dr 2 . . . dr N (8.10) 

can be inverted: given a groundstate density po( r ) h is possible, in principle, 
to calculate the corresponding ground-state wave function \I / o( r i ) r 2 --j y n)- 
This means that is a functional of po- Consequently, all ground-state 
observables are functionals of p , too. If \l/ can be calculated from p and 
vice versa, both functions are equivalent and contain exactly the same in- 
formation. At first sight this seems impossible: how can a function of one 
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(vectorial) variable r be equivalent to a function of N (vectorial) variables 
ri, r 2 ..., rjv? How can one arbitrary variable contain the same information as 
N arbitrary variables? 

The crucial fact which makes this possible is that knowledge of po(r) 
implies implicit knowledge of much more than that of an arbitrary func- 
tion /(r). The ground-state wave function \l/ must not only reproduce the 
ground-state density, but also minimize the energy. For a given ground-state 
density Po(r), we can write this requirement as 

E vfi =min{V\f + Vee + VeN\*) (8.11) 

where E Vt o denotes the ground-state energy in potential v(r). The preceding 
equation tells us that for a given density po(r) the ground-state wave function 
is that which reproduces this po(r) and minimizes the energy. For an 
arbitrary density p(r), we define the functional 

E v \p]=min{*\f + V ee + V eN \V). (8.12) 

w— >p 

If p is a density different from the ground-state density po i n potential v(r), 
then the \I/ that produce this density p are different from the ground-state 
wave function ty , and according to the variational principle the minimum 
obtained from E v [p] is higher than (or equal to) the ground-state energy 
E Vj0 = E v [po\. Thus, the functional E v [p] is minimized by the ground-state 
density po, and its value at the minimum is E v $. 
The total energy functional can be written as 

E v [p] = min(V\f + Vjp) + d 3 r p(r)v(r) (8.13) 
=:F\p] + V\p] (8.14) 

where the internal-energy functional F[p] = min{^\T -^Ve^) is independent 

of the potential v(r), and thus determined only by the structure of the oper- 
ators T and V ee . This universality of the internal-energy functional allows us 
to define the ground-state wave function as that antisymmetric iV-particle 
function that delivers the minimum of F[p] and reproduces p Q . If the ground 
state is nondegenerate, this double requirement uniquely determines $ m 
terms of po(r), without having to specify v(r) explicitly. 4 Equations 8.11 
to 8.14 constitute the constrained-search proof of the Hohenberg-Kohn the- 
orem, given independently by M. Levy and E. Lieb. The original proof by 
Hohenberg and Kohn proceeded by assuming that was not determined 
uniquely by po an d showed that this produced a contradiction to the varia- 



4 Note that this is exactly the opposite of the conventional prescription to specify the 
Hamiltonian via v(r), and obtain from solving Schrodinger equation, without having to 
specify p(r) explicitly. 
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tional principle. Both proofs, by constrained search and by contradiction, are 
elegant and simple. In fact, it is a bit surprising that it took 38 years from 
Schrodinger's first papers on quantum mechanics to Hohenberg and Kohn's 
1964 paper containing their famous theorem. 



8.3.1 Meaning and implications of the Hohenberg- Kohn 
theorem 

Here we provide a commented summary of the content of the Hohenberg- 
Kohn (HK) theorem. 

(1) The nondegenerate ground-state (GS) wave function is a unique func- 
tional of the ground state (GS) density 5 

*o(r,r 2j ...,r JV ) = *[p (r)]. (8.15) 

This is the essence of the HK theorem (often called the first Hohenberg- 
Kohn theorem). As a consequence, the GS expectation value of any ob- 
servable O is a functional of Po( r ) ; too 

Oo = 0[ Po ] = (^M|0|^M>. (8.16) 

(2) For the special case in which the observable is the energy of the system, 

E v , = O[p ] = E v [p ] = m Po }\H\^[ Po }) . (8.17) 

where, H = T + V ee + V e N , the following variational property holds 

EM <E v \p'} (8.18) 

where p is GS density in potential V and pf is some other density. This 
is very similar to the usual variational principle for wave functions. From 
a calculation of the expectation value of a Hamiltonian with a trial wave 
function \P' that is not its GS wave function \Po one can never obtain an 
energy below the true GS energy, 

E vfi = E V [V ] = (*o|£|*o> < (y'\H\y'} = E V [V] (8.19) 



5 If the ground state is degenerate, several of the degenerate ground-state wave functions 
may produce the same density, so that a unique functional ^f[p] does not exist, but by 
definition these wave functions all yield the same energy, so that the functional E v [p] 
continues to exist and to be minimized by p$. A universal functional F[p] can also still be 
defined. 
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Similarly, in exact DFT, if E[p] for fixed v ext is evaluated for a density that 
is not the GS density of the system in potential v ext , one never finds a result 
below the true GS energy. This is what 8.18 says, and it is so important 
for practical applications of DFT that it is sometimes called the second 
Hohenberg-Kohn theorem. 

In performing the minimization of E v [p] the constraint that the total par- 
ticle number N is an integer is taken into account by means of a Lagrange 
multiplier, replacing the constrained minimization of E v [p] by an uncon- 
strained one of E v [p] — pN. Since N = J d 3 r p(r), this leads to 

Sp(r) dN v ' 

where p is the chemical potential. 



(3) Recalling that the kinetic and interaction energies of a nonrelativistic 
Coulomb system are described by universal operators, we can also write E v 

as 

E v [p] = T[p) + V ee [p] + V eN [p] = F[p] + V eN [p] , (8.21) 

where T[p] and V ee [p] are universal functionals (defined as expectation values 
of the type of Eq. 8.17), independent of Ugjv(r) = v ex t(r) (sometimes simply 
called v(r)). On the other hand, the potential energy in a given potential 
v(r) is the expectation value of the potential 

^w = E^^ ( 8 - 22 ) 

which reads 

V ext [p] = J d 3 r p{v)v ext {r) (8.23) 

and is obviously nonuniversal (it depends on v ext (r), i.e., on the system under 
study), but very simple: once the system is specified, i.e., v ext (r) is known, 
the functional Kxt[p] is known explicitly. 



8.3.2 Density- Functional Theory in practice 

After these abstract considerations let us now consider one way in which one 
can make practical use of DFT. Assume we have specified our system (i.e., 
f(r) is known). Assume further that we have reliable approximations for 
V ee [p] and T[p\. In principle, all one has to do then is to minimize the sum 
of kinetic, interaction and potential energies 

E v [p] = T[p] + V ee [p] + V eN [p] = T[p] + V ee [p] + / d 3 r p{r)v ext (r) (8.24) 
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with respect to p(r). The minimizing function po( r ) is the system's GS 
charge density and the value E v< q = E v [p ] is the GS energy. Assume now 
that v ex t(r) depends on a parameter a. This can be, for example, the lattice 
constant in a solid or the angle between two atoms in a molecule. Calculation 
of E v> q for many values of a allows one to plot the curve E v ^ (a) and to find 
the value of a that minimizes it. This value, a , is the GS lattice constant or 
angle. In this way one can calculate quantities like molecular geometries and 
sizes, lattice constants, unit cell volumes, charge distributions, total energies, 
etc. 

By looking at the change of E Vj0 (a) with a one can, moreover, calculate 
compressibilities, phonon spectra and bulk moduli (in solids) and vibrational 
frequencies (in molecules). By comparing the total energy of a composite 
system (e.g., a molecule) with that of its constituent systems (e.g., individual 
atoms) one obtains dissociation energies. By calculating the total energy 
for systems with one electron more or less one obtains electron affinities 
and ionization energies. Using the Hellman-Feynman theorem 6 one can 
calculate forces on atoms from the derivative of the total energy with respect 
to the nuclear coordinates. All this follows from DFT without having to 
solve the many-body Schrodinger equation and without having to make a 
single-particle approximation. 

In theory it should be possible to calculate all observables, since the 
HK theorem guarantees that they are all functionals of po- In practice, 
one does not know how to do this explicitly. Another problem is that the 
minimization of E v [p\ is, in general, a tough numerical problem on its own. 
And, moreover, one needs reliable approximations for T[p] and V ee [p] to begin 
with. In the next section, on the Kohn-Sham equations, we will see one widely 
used method for solving these problems. Before looking at that, however, it 
is worthwhile to recall an older, but still occasionally useful, alternative: the 
Thomas-Fermi approximation. 



6 For the case of a = R/ (the coordinate of a nucleus I), the Hellmann-Feynman theorem 
says that, for the GS wavefunction or GS density, 

p d(E GS ) 

/f)ff 
d 3w vl/*( ri;r2 ,..., rjv )__vl/( ri , r2 ,..., rjv ) 

= G>K,o[p](Rj) 

<9R 7 
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8.3.3 The Thomas- Fermi approximation and the local 
density approximation (LDA) 

In this approximation one sets 

V-«V B = \jfrl&<!Mn (8.25) 

i.e., approximates the full interaction energy by the Hartree energy, the elec- 
trostatic interaction energy of the charge distribution p(r). One further ap- 
proximates, initially, 

T[p\ « T LDA [p] = J t hom {p{r))d 3 r (8.26) 

where t hom (p) is the kinetic-energy density of a homogeneous interacting 
system with (constant) density p. Since it refers to interacting electrons 
t hom (p) is not known explicitly, and further approximations are called for. 
As it stands, however, this formula is already a first example of a local 
density approximation (LDA). In this type of approximation one imagines 
the real inhomogeneous system (with density p(r) in potential v(r)) to be 
decomposed in small cells in each of which p(r) and v(r) are approximately 
constant. In each cell (i.e., locally) one can then use the per- volume energy 
of a homogeneous system to approximate the contribution of the cell to the 
real inhomogeneous one. Making the cells infinitesimally small and summing 
over all of them yields Eq. 8.26. 

For a noninteracting system (specified by subscript s, for 'single-particle') 
the function t^ om (p) is known explicitly, 

t h s om ( p ) = 3/l 2 (37r 2 ) 2/3 p 5/ 7(i0m e ) ( 8 .27) 

We can therefore use the following approximation 

T[p] « T LDA [p] » T s LDA [p] = J t h np(v))dh (8.28) 

where T^ DA [p] is the local-density approximation to T[p\. The Thomas-Fermi 
approximation then consists in combining Eq. 8.25 with Eq. 8.26 and writing 

E[p] = T[p] + V ee [p] + V eN [p] « E TF [p] = T s LDA [p] + V H [p] + V eN [p}. (8.29) 

A major defect of the Thomas- Fermi approximation is that within it molecules 
are unstable: the energy of a set of isolated atoms is lower than that of 
the bound molecule. This fundamental deficiency, and the lack of accu- 
racy resulting from neglect of correlations in 8.25 and from using the local 
approximation 8.28 for the kinetic energy, limit the practical use of the 
Thomas-Fermi approximation in its own right. However, it is found a most 
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useful starting point for a large body of work on improved approximations in 
chemistry and physics. More recent approximations for T[p] can be found in 
the context of orbital-free DFT. The extension of the local-density concept 
to the exchange-correlation energy is at the heart of many modern density 
functionals. 



8.4 Kohn-Sham Density Functional Theory 

Density-functional theory can be implemented in many ways. The mini- 
mization of an explicit energy functional, discussed up to this point, is not 
normally the most efficient among them. Much more widely used is the Kohn- 
Sham approach. Interestingly, this approach owes its success and popularity 
partly to the fact that it does not exclusively work in terms of the particle (or 
charge) density, but brings a special kind of wave functions (single-particle 
orbitals) back into the game. As a consequence, DFT then looks formally 
like a single-particle theory, although many-body effects are still included via 
the so-called exchange-correlation functional. We will now see in some 
detail how this is done. 



8.4.1 The exchange-correlation energy 

The Thomas-Fermi approximation for T[p] is not very good. A more accurate 
scheme for treating the kinetic-energy functional of interacting electrons, 
T[p], is based on decomposing T[p] into one part that represents the kinetic 
energy of noninteracting particles of density p, i.e., the quantity called above 
T s [p], and one that represents the remainder, denoted T c [p] (the sub-scripts 
s and c stand for 'single-particle' and 'correlation', respectively) 7 . 

T[p]=T s [p}+T c [p}. (8.30) 

T s [p] is not known exactly as a functional of p, but it is easily expressed in 
terms of the single-particle orbitals 4>i(r) of a noninteracting system with 
the density p, as 

N 

T M = ~Y1 J d'rtf (r)V%(r) (8.31) 

because for noninteracting particles the total kinetic energy is just the sum 
of the individual kinetic energies! 



7 T S is defined as the expectation value of the kinetic-energy operator T with the Slater 
determinant (\p) defined by the density p, i.e., T s [p] = (G>[p]\T\&[p\) . Similarly, the full 
kinetic energy is defined as T[n] = (^[p]\T\^[p]) . All consequences of antisymmetrization 
(i.e., exchange) are described by employing a determinantal wave function in defining T s . 
Hence, T c , the difference between T s and T is a pure correlation effect. 
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Since all <f>i(r) are functionals of p (p(r) = l<M r )| 2 ); this expression 
for T s is an explicit orbital functional but an implicit density functional, 
T s [p] = T s [<f)i[p}}, where the notation indicates that T s depends on the full set 
of occupied orbitals 0j, each of which is a functional of p. We now rewrite 
the exact energy functional as 

E[p] = T[p] + V ee [p] + V eN [p] = TsiMp]}] + V H [p] + EM + V e M, (8-32) 

where by definition E xc contains the differences T — T s (i.e. T c ) and V ee — Vh- 
This definition shows that a part of the correlation energy E c is due to the 
difference T c between the noninteracting and the interacting kinetic energies. 
Unlike the Thomas-Fermi equation, Eq. 8.36 is formally exact, but of course 
the functional E xc is unknown. 

The functional E xc [p] is called the exchange-correlation (xc) energy. 

It is often decomposed into E xc — E x + E c , where E x is due to the Pauli 
principle (exchange energy) and E c is due to electron correlation. (T c is then 
a part of E c .) The exchange energy can be written explicitly in terms of the 
single-particle orbitals as 

EJLitJA}] = ~E /<*/ (8.33) 



8.4.2 The Kohn-Sham equations 

Since T s is now written as an orbital functional one cannot directly minimize 
Eq. 8.36 with respect to p. Instead, one commonly employs a scheme sug- 
gested by Kohn and Sham for performing the minimization indirectly. 
This scheme starts by writing the minimization as 

= (8.34) 

ST s [p] 5V ext [ P ] 5V H [p] 5EM / R «x 
6p(r) 5p(r) Sp(r) 6p(r) [ } 

= ^4 + v ^ r ) + + v ^ ( 8 ' 36 ) 
5p(r) 

As a consequence of Eq. 8.22, 5V ext /Sp(r) = v ext (r), the 'external' potential 
the electrons move in (we use V ext for V c n). The term 5Vh/p{y) simply yields 
the Hartree potential. For the term 5E xc /5p, which can only be calculated 
explicitly once an approximation for E xc has been chosen, one commonly 
writes v xc . 

Consider now a system of noninteracting particles moving in the potential 
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v s (r). For this system the minimization condition is simply 



= (8.37) 
5p(r) 

- S ™ + S M (8.38) 



Sp(r) Sp(r) 
5p(r) 



(8.39) 



since there are no Hartree and xc terms in the absence of interactions. The 
density solving this Euler equation is p s (r). Comparing this with Eq. 8.36 
we find that both minimizations have the same solution p s (r) = p(r), if v s is 
chosen to be 

v s (r) = v ext (r) + v H (r) + v xc (r) . (8.40) 

Consequently, one can calculate the density of the interacting (many-body) 
system in potential v ext (r), described by a many-body Schrodinger equation 
of the form 8.1, by solving the equations of a noninteracting (single-body) 
system in potential v s (r). 8 

In particular, the Schrodinger equation of this auxiliary system, 



v 2 

-y + 



eMv) (8.4i; 



yields orbitals that reproduce the density p(r) of the original system (these 
are the same orbitals employed in Eq. 8.31), 



N 



p(r)^p s (r) = ^/# ( r)| 2 , (8.42) 

i 

where /, is the occupation of the i th orbital. 9 Eqs. 8.40 to 8.42 are the cel- 
ebrated Kohn-Sham (KS) equations. They replace the problem of mini- 
mizing E[p] by that of solving the Schrodinger equation for a non-interacting 
system. 



8 The question whether such a potential v s (r) always exists in the mathematical sense 
is called the noninteracting u-representability problem. It is known that every interacting 
ensemble u-representable density is also noninteracting ensemble v-representable, but only 
in discretized systems has it been proven that all densities are interacting ensemble v- 
representable. It is not known if interacting ensemble-representable densities may be 
noninteracting pure-state representable (i.e, by a single determinant), which would be 
convenient (but is not necessary) for Kohn-Sham calculations. 

9 Normally, the occupation numbers /j follows an Aufbau principle (Fermi statistics) 
with fi = 1 for i < N, /, = for i > N, and < fi < 1 for i = N (i.e., at most the 
highest occupied orbital can have fractional occupation). 
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Since both i>h and v xc depend on p, which depends on the fa, which in 
turn depend on v s , the problem of solving the KS equations is a nonlinear 
one. The usual way of solving such problems is to start with an initial guess 
for p(r), calculate the corresponding v s (r), and then solve the differential 
equation 8.41 for the fa. From these one calculates a new density, using 8.42, 
and starts again. 

The process is repeated until it converges. The technical name for this 
procedure is self-consistency cycle, which is the same we have encountered 
in the solution of the HF SC field equations. Different convergence criteria 
(such as convergence in the energy, the density, or some observable calcu- 
lated from these) and various convergence-accelerating algorithms (such as 
mixing of old and new effective potentials) are in common use. Only rarely 
it requires more than a few dozen iterations to achieve convergence, and even 
rarer are cases where convergence seems unattainable, i.e., a self-consistent 
solution of the KS equation cannot be found. 

Once one has a converged solution po, one can calculate the total energy 
from Eq. 8.36 or, equivalently and more conveniently, from 

£o = f>~ fd 3 r I d 3 r / PoWMp _ I d 3 rVxc{r)po{r)+Exc[po] (8 . 43) 
i=1 J J I I J 

Equation 8.43 follows from writing V ext [p] in 8.36 by means of 8.40 as 

V ext [p\ = J d 3 v ext (r)p(r) 

= J d 3 r [v s (r) -v H {r) - v xc (r)]p(r) 

= V r [p] - J d 3 r[v H (r) + v xc (r)] 

and identifying the energy of the noninteracting (Kohn-Sham) system as 
E s = J2? =1 e i = T s + V s . 

The eigenvalues of the Kohn-Sham equations 

Equation 8.43 shows that E is not simply the sum of all In fact, it 
should be clear from our derivation of Eq. 8.41 that the are introduced 
as completely artificial objects: they are the eigenvalues of an auxiliary one- 
body equation whose eigenfunctions (orbitals) yield the correct density. It is 
only this density that has strict physical meaning in the KS equations. The 
KS eigenvalues, on the other hand, in general bear only a semiquantitative 
resemblance with the true energy spectrum, but are not to be trusted quan- 
titatively. 
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The main exception to this rule is the highest occupied KS eigenvalue. 
Denoting by En(M) the N th eigenvalue of a system with M electrons, one 
can show rigorously that En{M) = —I, the negative of the first ionization 
energy of the N-body system, and En+i{N + 1) = —A, the negative of the 
electron affinity of the same iV-body system. These relations hold for the 
exact functional only. When calculated with an approximate functional of 
the LDA or GGA type (see below), the highest eigenvalues usually do not 
provide good approximations to the experimental I and A. Better results 
for these observables are obtained by calculating them as total-energy differ- 
ences, according to I = E (N - 1) - E (N) and A = E (N) - E (N + 1), 
where E Q (N) is the ground-state energy of the iV-body system. Alternatively, 
self-interaction corrections (SIC) can be used to obtain improved ionization 
energies and electron affinities from Kohn-Sham eigenvalues. 
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Figure 8.1: Schematic description of some important Kohn-Sham 
eigenvalues relative to the vacuum level, denoted by 0, and their 
relation to observables. See main text for explanations. 



Figure 8.1 illustrates the role played by the highest occupied and lowest 
unoccupied KS eigenvalues, and their relation to observables. For molecules, 
HOMO(A r ) is the highest-occupied molecular orbital of the iV-electron sys- 
tem, HOMO(A r + 1) that of the N + 1-electron system, and LUMO(A r ) the 
lowest unoccupied orbital of the iV-electron system. In solids with a gap, 
the HOMO and LUMO become the top of the valence band and the bot- 
tom of the conduction band, respectively, whereas in metals they are both 
identical to the Fermi level. The vertical lines indicate the Kohn-Sham (sin- 
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gleparticle) gap A#s, the fundamental (many-body) gap A, the derivative 
discontinuity of the xc functional, A xc , the ionization energy of the interact- 
ing iV-electron system I(N) = —En(N) (which is also the ionization energy 
of the Kohn-Sham system Iks(N)), the electron affinity of the interacting 
Nelectron system A(N) = —en+i(N + 1) and the Kohn-Sham electron affin- 
ity A KS (N) = -e N+1 (N). 

Given the auxiliary nature of the other Kohn-Sham eigenvalues, it comes 
as a great (and welcome) surprise that in many situations (typically charac- 
terized by the absence of strong correlations) the Kohn-Sham eigenvalues 
do, empirically, provide a reasonable first approximation to the actual energy 
levels of extended systems and molecules. 



8.5 Making DFT practical: Approximations 

The basic approximation in DFT consists in the construction of an expres- 
sion for the unknown xc functional E xc [p], which contains all many-body 
aspects of the problem. This chapter is intended to give you an idea of 
what types of functionals exist, and to describe what their main features are, 
separately for local functionals (TF, LDA and X Q ), semilocal, or gradient- 
dependent functionals (GEA and GGA), and nonlocal functionals (hybrids, 
orbital functionals such as meta-GGAs, EXX and self interaction corrected 
SIC ones. This chapter does deal only most superficially with the actual con- 
struction of these functionals. For more details on functional construction 
and testing the students are referred to the reviews on this subject (see for 
instance the literature given in chapter 8.1). 



8.5.1 Local functionals: LDA 

Historically (and in many applications also practically) the most important 
type of approximation is the local-density approximation (LDA). To under- 
stand the concept of an LDA recall first how the noninteracting kinetic en- 
ergy T s [p] is treated in the Thomas- Fermi approximation: In a homogeneous 
system one knows that, per volume 

t h s ° m (p) = To^( 37r2 ) 2/3 P 5/3 ( 8 - 44 ) 

where p = const. In an inhomogeneous system, with p = p(r), one approxi- 
mates locally 

*.(r) « tT(p(r)) = -^-(3tt 2 ) 2 / 3 p(r) 5 / 3 (8.45) 
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and obtains the full kinetic energy by integration over all space 



T S LDA [ P ] = J d*rt h r(p(r)) = ^(Svr 2 ) 2 / 3 J d 3 r p(rf* (8.46) 

For the kinetic energy the approximation T s [p] m Tf LDA [p] is much inferior 
to the exact treatment of T s in terms of orbitals, offered by the Kohn-Sham 
equations, but the LDA concept turned out to be highly useful for another 
component of the total energy, the exchange-correlation energy E xc [p\. For 
the exchange energy E x [p] the procedure is very simple, since the per volume 
exchange energy of the homogeneous electron liquid is known exactly 

D hom(\ _ 3ql 1 3 \ V 3 4/3 



0) = -j(-J P" A (8-47) 



so that 

-,2 



E" DA lp] = - 3j t(l) 1/3 Jd 3 rp 4/3 (8-48) 



This is the LDA for E T . 



For the correlation energy E c [p] the situation is more complicated since 
e c° m (p) is n °t known exactly: the determination of the correlation energy of 
a homogeneous interacting electron system (an electron liquid) is already a 
difficult many-body problem on its own! Early approximate expressions for 
e c° m (p) were based on applying perturbation theory (e.g. the random phase 
approximation) to this problem. These approximations became outdated 
with the advent of highly precise Quantum Monte Carlo (QMC) calculations 
for the electron liquid, by Ceperley and Alder (D. M. Ceperley and B. J. 
Alder, Phys. Rev. Lett., 45, 566 (1980)). Modern expressions for e^ om (p) 
are parametrizations of these data. These expressions are implemented in 
most standard DFT program packages and in typical applications give al- 
most identical results. On the other hand, the earlier parametrizations of 
the LDA, based on perturbation theory, can occasionally deviate substan- 
tially from the QMC ones, and are better avoided. 

Independently of the parametrization, the LDA for E xc consists in 
E xc [p] « E™ A \p] = J e h x r(p)\^ P (r) d 3 r = J eJr(p(r)) d 3 r (8.49) 
where e x ° m = e x om + e h ° m . The corresponding xc potential 



LDA 



MM 



de h x r(p) 



dp 



(8.50) 



This approximation for E xc [p] has proved amazingly successful, even when 
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applied to systems that are quite different from the electron liquid that forms 
the reference system for the LDA. A partial explanation for this success of 
the LDA is systematic error cancellation: typically, LDA underestimates E c 
but overestimates E x , resulting in unexpectedly good values of E xc . 

For many decades the LDA has been applied in, e.g., calculations of band 
structures and total energies in solid-state physics. In quantum chemistry it is 
much less popular, because it fails to provide results that are accurate enough 
to permit a quantitative discussion of the chemical bond in molecules (so- 
called 'chemical accuracy' requires calculations with an error of not more 
than about lkcal/mol = 0.04336eV/ particle). At this stage it may be worth- 
while to recapitulate what practical DFT does, and where the LDA enters its 
conceptual structure: What real systems, such as atoms, molecules, clusters 
and solids, have in common, is that they are simultaneously inhomogeneous 
(the electrons are exposed to spatially varying electric fields produced by the 
nuclei) and interacting (the electrons interact via the Coulomb interaction). 
The way density-functional theory, in the local-density approximation, deals 
with this inhomogeneous many-body problem is by decomposing it into two 
simpler (but still highly nontrivial) problems: the solution of a spatially ho- 
mogeneous interacting problem (the homogeneous electron liquid) yields the 
uniform xc energy e^° m (p), and the solution of a spatially inhomogeneous 
noninteracting problem (the inhomogeneous electron gas described by the 
KS equations) yields the particle density. Both steps are connected by the 
local-density potential ( 8.50), which shows how the xc energy of the uniform 
interacting system enters the equations for the inhomogeneous noninteract- 
ing system. 



8.5.2 Semilocal functionals: GEA, GGA and beyond 

In the LDA one exploits knowledge of the density at point r. Any real system 
is spatially inhomogeneous, i.e., it has a spatially varying density p(r), and it 
would clearly be useful to also include information on the rate of this variation 
in the functional. A first attempt at doing this were the so-called gradient- 
expansion approximations (GEA). In this class of approximation one tries 
to systematically calculate gradient-corrections of the form | Vp(r)|, | Vp(r)| 2 , 
V 2 p(r), etc., to the LDA. A famous example is the lowest-order gradient 
correction to the Thomas- Fermi approximation for T s [p], 

Ts[p] - Tf [ P ] = Tt DA [p] + ^Jd*r (8.51) 
This second term on the right-hand side is called the Weizsacker term. 
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Similarly, in 



E x [p] « E^[p] = E™ A [p] 



We 




V 



|Vp(r)| 2 
p(r)4/3 



(8.52) 



432 7 r(37r 2 ) 1 /3 



the second term on the right-hand side is the lowest-order gradient correction 
to E^ DA . In practice, the inclusion of low-order gradient corrections almost 
never improves on the LDA, and often even worsens it. Higher-order correc- 
tions (e.g., oc |Vp(r)| Q or V 6eta p(r) with a,/3 > 2), on the other hand, are 
exceedingly difficult to calculate, and little is known about them. 



In this situation it was a major breakthrough when it was realized, in the 
early eighties, that instead of power-series-like systematic gradient expan- 
sions one could experiment with more general functions of p(r) and Vp(r), 
which need not proceed order by order. Such functionals, of the general form 



have become known as generalized-gradient approximations (GGAs). 
Different GGAs differ in the choice of the function f(p, Vp). Note that this 
makes different GGAs much more different from each other than the different 
parametrizations of the LDA: essentially there is only one correct expression 
for e£° m (p), and the various parametrizations of the LDA are merely different 
ways of writing it. On the other hand, depending on the method of construc- 
tion employed for obtaining /(p, Vp) one can obtain very different GGAs. 
In particular, GGAs used in quantum chemistry typically proceed by fitting 
parameters to test sets of selected molecules. On the other hand, GGAs used 
in physics tend to emphasize exact constraints. Nowadays the most popu- 
lar (and most reliable) GGAs are PBE (denoting the functional proposed in 
1996 by Perdew, Burke and Ernzerhof) in physics, and BLYP (denoting the 
combination of Becke's 1988 exchange functional with the 1988 correlation 
functional of Lee, Yang and Parr) in chemistry. Many other GGA-type func- 
tionals are also available, and new ones continue to appear. 




(8.53) 
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method 


-E [a.u. 


Thomas-Fermi 


625.7 


Hartree-Fock 


526.818 


LDA (exchange only) 


524.517 


LDA (exch. and corr.) 


525.94 


LDA-SIC(PZ) 


528.393 


GGA BLYP 


527.551 


experiment 


527.6 



Table 8.1: Ground state energy in atomic units (1 a.u.= 1 Hartree = 2 
Rydberg = 27.21 eV = 627.5 kcal/mol) of the Ar atom (Z=18), obtained 
with HF and DFT for some representative functionals. 

Quite generally, current GGAs seem to give reliable results for all main 
types of chemical bonds (covalent, ionic, metallic and hydrogen bonds). For 
van der Waals interactions, however, common GGAs and LDA fail. To de- 
scribe these very weak interactions several more specialized approaches have 
been developed within DFT. Both in physics and in chemistry the widespread 
use of GGAs has lead to major improvements as compared to LDA. Chem- 
ical accuracy, as defined above, has not yet been attained, but is not too 
far away either. A useful collection of explicit expressions for some GGAs 
can be found in the literature (see for instance P. Ziesche, S. Kurth and J. P. 
perdew, Comp. Mat. Sc% 11, 122 (1998)). 

No systematic attempt at comparing explicit functionals can be made here, 
but many detailed comparisons are available in the literature. For purely 
illustrative purposes only, Table 8.1 contains ground-state energies of the 
Ar atom, obtained with several of the methods discussed in this and in the 
following sections. 

8.5.3 Orbital functionals and other nonlocal approxima- 
tions: hybrids, Meta-GGA, SIC, etc. 

In spite of these advances, the quest for more accurate functionals goes ever 
on, and both in chemistry and physics various beyond-GGA functionals have 
appeared. Perhaps the most popular functional in quantum chemistry is 
B3LYP. This is a combination of the LYP GGA for correlation with Becke's 
three-parameter hybrid functional B3 for exchange. Common hybrid func- 
tionals, such as B3, mix a fraction of Hartree-Fock exchange 

E41MP]}] =-\Tj *>7 d V ^ r) ^y* (r,) (8.54) 

into the DFT exchange functional (other mixtures are also possible). The 
construction of hybrid functional involves a certain amount of empiricism in 
the choice of functionals that are mixed and in the optimization of the weight 



105 



factors given to the HF and DFT terms. Formally, this might be considered a 
drawback, but in practice B3 has proven to be the most successful exchange 
functional for chemical applications, in particular when combined with the 
LYP GGA functional for E c . A more extreme example of this semiempirical 
mode of construction of functionals is the Becke's 1997 hybrid functional, 
which contains 10 adjustable parameters. 

Another recent beyond-GGA development is the emergence of so-called Meta- 
GGAs, which depend, in addition to the density and its derivatives, also on 
the Kohn-Sham kinetic-energy density r(r) 



so that E xc can be written as E xc [p(r), Vp(r), r(r)]. The additional degree of 
freedom provided by r is used to satisfy additional constraints on E xc , such 
as a self-interaction-corrected correlation functional and a finite exchange 
potential at the nucleus. In several recent tests Meta-GGAs have given fa- 
vorable results, even when compared to the best GGAs, but the full potential 
of this type of approximation is only beginning to be explored systematically. 

As we have seen in the case of T s , it can be much easier to represent a func- 
tional in terms of single-particle orbitals (like the Kohn-Sham orbitals) than 
directly in terms of the density. Such functionals are known as orbital func- 
tionals, like for instance T s in Eq. 8.31. Another important orbital-dependent 
functional is the exchange energy (Fock term) of Eq. 8.54. The Meta-GGAs 
and hybrid functionals mentioned above are also orbital functionals, because 
they depend on the kinetic energy density (8.55), and on a combination of 
the orbital functional (Eq. 8.54) with ordinary GGAs, respectively. Still an- 
other type of orbital functional is the self-interaction correction (SIC). 
Most implementations of SIC make use of the following expression proposed 
by J. P. Perdew and A. Zunger (PZ-SIC), 



Eg GA ~ SIC ip a , p p ] = Eg GA [ Pa , p p ] - (EhIpu,] - Eg GA [p ia , 0]) (8.56) 



which subtracts, orbital by orbital, the contribution the Hartree and the xc 
functionals would make if there was only one electron in the system. Here 
a G {a,/3}. This correction can be applied on top of any approximate den- 
sity functional, and ensures that the resulting corrected functional satisfies 

Ex GA ~ SIC [p^ , 0] = —Eh[p^] for a one-electron system with density p^\ 
The LDA is exact for a completely uniform system, and thus is self-interaction 
free in this limit, but neither it nor common GGAs satisfy the requirement 
of freedom from self-interaction in general, and even Meta-GGAs have a re- 
maining self-interaction error in their exchange part. This self-interaction is 
particularly critical for localized states, such as the d states in transition- 
metal oxides. For such systems PZ-SIC has been shown to greatly improve 
the uncorrected LDA and GGA functionals, but for thermochemistry PZ-SIC 




(8.55) 
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does not seem to be significant. 

Unfortunately the PZ-SIC approach, which minimizes the corrected energy 
functional with respect to the orbitals, does not lead to Kohn-Sham equations 
of the usual form, because the resulting effective potential is different for each 
orbital. As a consequence, various specialized algorithms for minimizing the 
PZ-SIC energy functional have been developed. For finite systems, PZ-SIC 
has also been implemented by means of the "optimized effective potential" 
OEP method, which produces a common local potential for all orbitals. The 
study of this theory is however beyond the scope of this lecture. The inter- 
ested students will find more about this subject and other more advanced 
topics in the master course. 
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Chapter 9 



Ab-initio Molecular Dynamics in 
the Ground State 



Since an electron is much lighter than a nucleus, one may usually assume 
the Born-Oppenheimer (BO) adiabatic approximation. That is, one ignores 
the coupling between the motions of electrons and of nuclei and performs 
the integrations of these two motions separately at each time step. When 
the motion of nuclei is not very fast, and when the temperature is very 
low, one assumes also that the electronic states are at the ground state. 
The term, "adiabatic potential surface" or "Born-Oppenheimer surface" (BO 
surface) is often used to indicate an energy surface, in the 3iV-dimensional 
(atomic) configuration space (N refers to the number of atoms in the system), 
calculated using the BO approximation. 

Here we should make some remarks on the use of the BO approximation. 
First, for systems composed of light atoms such as hydrogen, there are cases in 
which, even in the adiabatic approximation, the zero-point vibration energy 
of the nuclei may become important in a discussion of the subtle energy 
differences between stable atomic configurations. In such a case, one has 
to include the contribution from the zero-point energy. Second, the BO 
approximation itself may break down if, for example, the electrons are not 
in the ground state or the atoms move at very high speeds. 

The resulting Schrodinger equation for the ionic motion is given by Eq. ?? 



^(R)=J5^(R) 



where y?"(R) is the amplitude for the nuclei to have positions R when the 
electrons are in the state ^ e n (R = (Ri, R2, •••))• 

The total electronic energy, e n (R), which includes the interaction among the 
nuclei, Vjvtv, plays the role of an effective potential energy for the nuclei. 



1 The index n labels the different solutions of the electronic Scrodinger equation, n = 
stays for the ground state. 
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In the quantum classical approximation only the electrons are quantized, 
while the nuclei are treated at a classical level (point charges). As a conse- 
quence the above Schrodinger equation for the nuclear dynamics is replaced 
by its classical equivalent, the Newton equation 

F I = m I R I , (9.1) 
where the forces, Fj, are computed from the total electronic potential, £ n (R) 

^n(R) 



<9R, 



■V/e„(R) (9.2) 



On the basis of density functional theory, the right-hand side of this 
equation can be written as 

-V 7 e n (R) = -V/ \ R Zl -n i " / pn ^ ^M\ T ~ R/ D dr ( 9 - 3 ) 
j&i) ' 1 A J 

where the first and the second terms represent, respectively, the Coulomb re- 
pulsive force between the nuclei and the Coulomb attractive force exerted on 
the nuclei by the electron cloud. Here, Vjir) denotes either a pseudopoten- 
tial, or —Zijr in a all electron approach (r = |r|), and p n (r) is the electronic 
density of the ground state (if n = 0) or a selected excited state (n > 1). 



9.1 The Hellmann-Feynman Forces 

An important ingredient in all dynamics methods is the efficient calculation 
of the forces acting on the nuclei. The straightforward numerical evaluation 
of the derivative 

Fj = -V/(tf |ftel|tto> ( 9 ' 4 ) 

in terms of finite-difference approximation of the total electronic energy is 
both too costly and too inaccurate for dynamical simulations. What happens 
if the gradients are evaluated analytically? 
In addition to the derivative of the Hamiltonian itself 

V / (*o|Kl|*0> = (*o|V/Wei|* ) + (V^olKil^o) + (tfolKilV/tfo) (9.5) 

there are in general also contributions from the variations of the wavefunc- 
tions (V/^o)- in general means here that these contributions vanish exactly 
if the wavefunction is an exact eigenfunction (or stationary state wavefunc- 
tion) of the Hamiltonian under consideration and the nuclear forces become 

Ff F = -(tfolV/fteil^o) (9.6) 
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This is the content of the often cited Hellmann-Feynman Theorem, which is 
also valid for many variational wavefunctions (e.g. Hartree-Fock or Kohn- 
Sham wavefunctions) . 



Ill 



9.2 Ehrenfest Molecular Dynamics 



The more straightforward molecular dynamical scheme amounts to compute 
the Ehrenfest force by actually solving numerically the coupled set of nuclear 
and electronic propagation equations: 



M 7 R 7 (t) 



ih- 



dt 



-V/(tf |ftel|*0> 

h 2 



(9.7) 
(9.8) 



The Ehrenfest approach describes a real time-dependent evolution of the 
electronic degrees of freedom and therefore includes rigorously non-adiabatic 
effects like transitions between different electronic states and 2 , which 
however are not considered in this chapter. 

Ehrenfest dynamics is just introduced here as a rigorous way of propagat- 
ing both nuclear and electronic degrees of freedom in the adiabatic regime. 
The application of the Ehrenfest molecular dynamics approach in compu- 
tational chemistry is however limited by the size of the propagation time 
step, which is of the order of 0.1 a.u., and therefore 10 to 100 smaller than 
the one used in Car-Parrinello and Born-Oppenheimer molecular dynamics, 
respectively. 



9.3 Born-Oppenheimer Molecular Dynamics 

The simplest way of including electronic structure in molecular dynamics 
simulations consists in straightforwardly solving the static electronic struc- 
ture problem in each molecular dynamics step given the set of fixed nuclear 
positions at the instant of time. Thus, the electronic structure part is re- 
duced to solving the time-independent quantum problem, e. g., by solving 
the time-independent Schodinger equation, concurrently to propagating the 
nuclei via classical molecular dynamics. In this picture, the time-dependence 
of the electronic structure is a consequence of the nuclear motion, and not 
explicit as for instance in Ehrenfest dynamics. 

The resulting Born-Oppenheimer molecular dynamics method is defined 

by 

MjKjit) = -V / min{(^ |^ei|*o)} (9.9) 
E V = H e ^ (9.10) 

for the electronic ground state. A deep difference with respect to Ehrenfest 
dynamics concerning the nuclear equation of motion is that the minimum of 



2 see for instance: J. C. Tully, in Modern Methods for Multidimensional Dynamics 
Computations in Chemistry, ed. D.L. Thompson (World Scientific, Singapore, 1998). 
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(T~i e i) has to be reached in each Born-Oppenheimer molecular dynamics step. 
In Ehrenfest dynamics, on the other hand, a wavefunction that minimized 
(Hei) initially will also stay in its respective minimum as the nuclei move. 



9.4 Car-Parrinello Method 

Since the paper by Car and Parrinello 3 first appeared in 1985, ab initio MD 
techniques have developed very rapidly. In their paper, Car and Parrinello 
proposed to integrate a Newtonian equation of motion for both the wavefunc- 
tions and the atomic coordinates. Such an approach was a completely new 
idea, beyond the previous common-sense idea that the electronic structure 
must be calculated by means of matrix diagonalization and self-consistent 
iterations. 

Since the basic equation of motion in the Car-Parrinello method is derived 
from a Lagrangian, total energy conservation is automatically guaranteed in 
the formulation (the total energy in this method can be considered as the true 
total energy of the system if ji, defined below, is sufficiently small), and their 
approach is especially suitable for microcanonical simulations (with constant 
total energy). Their subsequent very active studies with various applications 
demonstrated the validity and effectiveness of their approach. With this 
technique, it is no longer necessary to treat huge eigenvalue problems and it 
becomes possible to reduce significantly both the computational memory and 
the computational time. The reason for the recent very rapid development 
of ab initio MD techniques is mainly due to this approach. 

The basic equation of the Car-Parrinello method is given by the La- 
grangian £ as a functional of the wavefunction (=Kohn-Sham orbitals) ipi. 
(I denotes the level number of the electronic states) and the atomic position 
Rr (7 stands for the atom index) 

£ = Ef [ d 3 r\^\ 2 + 1 -J2 M i\ ji i\ 2 - e[m,{Ri}] (9.H) 
i J i 

Here, \i is a fictitious electron mass which governs the motion of the wave- 
functions, and e [{^i}, {R-/}] includes the electronic potential and the classi- 
cal interaction between the nuclei. Note that this Lagrangian coincides with 
the true Lagrangian of the system only in the limit fi — > 0, although we set 
/i finite. The important meaning of the existence of a Lagrangian is that 
the system has no energy dissipation and the total energy, including the fic- 
titious kinetic energy (KE) Hi\ijji\ 2 of the electrons, is conserved. Therefore, 
the atomic oscillation of this system is guaranteed to be preserved perma- 
nently without damping. From this Lagrangian, one may derive the equation 



3 R.Car and M. Parrinello, Phys. Rev. Lett, 55, 2493 (1985) 
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of motion for the wavef unctions 



fi^i = -n^ + J2 A ^u , (9.12) 

V 

where the K\ v are Lagrange multipliers to ensure the orthogonality of the 
wavefunctions ipi. This multiplier is determined so as to orthonormalize the 
wavef unctions. 

To ensure the Born-Oppenheimer adiabatic approximation in the Car-Parrinello 
formalism, one has to choose the left-hand side of Eq. (9.12) to be small 
enough. That is, the fictitious electron mass fi should be small. Here, one 
should comment on the fact that even with this choice of /i, the electronic 
wavefunctions oscillate very rapidly around the BO surface, following the 
nuclear motion. This rapid oscillation of the electronic states is of course 
not realistic, but it is a consequence of the second-order differentiation with 
respect to time in Eq. (9.12), like in the case of a wave equation. However, 
for sufficiently small /i, the amplitude of this oscillation is usually very small. 
On the other hand, the nuclear motion obeys the Newtonian equation 

d 2 

Mt—R^-VtE (9.13) 

with E being the sum of the Coulomb potential between the nuclei and the 
total energy of the electron system for nuclear positions, i.e. the physical 

total energy + Xl/(/ i /2)|^| 2 5 which is conserved. If the force on the right 
hand side of Eq. 9.13 is obtained, this equation can be integrated by means 
of, for example, the Verlet method. 



9.4.1 Why does the Car-Parrinello Method Works? 



In order to shed light on the title question, the dynamics generated by 
the Car-Parrinello Lagrangian Eq. 9.11 is analyzed in more detail invok- 
ing a "classical dynamics perspective" of a simple model system (eight sili- 
con atoms forming a periodic diamond lattice, local density approximation 
to density functional theory, normconserving pseudopotentials for core elec- 
trons, plane wave basis for valence orbitals, 0.3 fs time step with // = 300 
a.u., in total 20000 time steps or 6.3 ps; see G. Pastore et al., Phys. Rev. A 
44, 6334 (1991) for full details). 
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Figure 9.1: Vibrational density of states Eq. 9.14 (con- 
tinuous spectrum in upper part) and harmonic approxi- 
mation thereof (stick spectrum in lower part) of the elec- 
tronic degrees of freedom compared to the highest fre- 
quency phonon mode u;™ ax (triangle) for a model system. 



For this system the vibrational density of states or power spectrum of the 
electronic degrees of freedom, i.e. the Fourier transform of the statistically 
averaged velocity autocorrelation function of the classical fields 

POO 

f(cu) = / dtcos(ut) V(^(r,t)|^(r,0)> (9.14) 
Jo 

is compared to the highest frequency phonon mode a;™ ax of the nuclear sub- 
system in the nex figure. 

From this figure it is evident that for the chosen parameters the nuclear 
and electronic subsystems are dynamically separated: their power spectra 
do not overlap so that energy transfer from the hot to the cold subsystem is 
expected to be prohibitively slow. 
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Chapter 10 



Hybrid Quantum- 
Mechanics /Molecular-Mechanics 
Approach 

Quantum chemical methods are generally applicable and allow the calcula- 
tion of ground and excited state properties (molecular energies and struc- 
tures, energies and structures of transition states, atomic charges, reaction 
pathways etc.) 

Molecular Mechanical methods are restricted to the classes of molecules 
they have been designed for and their success strongly depends on the careful 
calibration of a large number of parameters. In addition chemical reactions 
cannot be described using Molecular Mechanics because the bonds are de- 
fined once and forever by the chosen force field. 

The development of the hybrid QM/MM approaches is guided by the 
general idea that large chemical systems may be partitioned into an elec- 
tronically important region which requires a quantum chemical treatment 
and a remainder which only acts in a perturbative fashion and thus admits 
a classical description. 

Using this approach it is possible to reduce the size of the quantum system, 
and treat the rest of the system, the environment, at a classical Molecular 
Mechanical level. In this way we can reduce substantially the amount of com- 
putational time and memory storage capacity required for the computation 
of complex systems made of thousands of atoms, like for instance a protein. 

The original molecular Hamiltonian for the QM system in the Born- 
Oppenheimer approximation 

tjQM 1 \ — v ^ t-72 1 v "v _2 V — "> Z[Z j ^ — \ Z\ ^ — \ 1 

tf« = "2 2^ Wj 1 ' 2 ^ Vn+ ^ |R 7 - R/| |Rj - r B | + ^ 

I 1 n I<J 11 Jl In 1 1 n| n<m 1 m 

(10.1) 
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is modified into 



" =-^m, ^■-2\.^ + l i w^7\-\ ! w^:\ + h 

-ErR^H+ETR^in (102) 

Kn KL 

where the indices K and L run over all MM atoms with charges Q. The last 
two terms in Eq. 10.2 describe the interaction of the electrons of the QM sub- 
system with the charges of the MM atoms, and the electrostatic interaction 
among the MM atoms. What is still missing - but can easily be added - 
is the part of the MM Hamiltonian dealing with the the molecular bonded- 
interactions and that you already know from classical Molecular Dynamics. 




Figure 10.1: Cartoon representing a simple QM/MM 
setup. The central QM region can become the active site 
of an enzyme, the MM region the protein matrix, and the 
boundary region the solvent. 
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